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First Year Report 
Kittisak Nui Tiyapan tf 
Introduction 


Geometry and mathematics 


The simplest geometrical figure is a circle while the simplest of 
all polygons is a triangle. The degree of freedom of triangles increases 
from the equilateral to the isosceles and the right triangles to the scalene 
triangles. While spending the summer of 1990 in a traineeship through 
AIESEC I was introduced to a geometrical puzzle which, as I came 
to learn later, is called the flexatube, but which I conjectured at the 
time that was from ancient China. This puzzle is made up of sixteen 
right isosceles triangles tiled into four squares, each comprising of four 
triangles, which are in turn joined together to form a loop. It can be 
easily made up using some hard papers, a pair of scissors and cello 
tape. There are in total twenty hinges, four of which are as long as 
the hypotenuse while the other sixteen have their length equal to the 
shorter side of the triangle. By turning these rigid triangles upon their 
hinges the inside surface of the strip can become the outside and vice 
versa. I found one solution a week later and back in Thailand a friend 
of mine found another solution. I think that these two are the only 
possible solutions but have not been able to proof it. 

Having this interest in geometrical puzzles I was delighted at the 
time to find that the symbol for the men’s toilets in Budapest is an 
equilateral triangle, while that for women’s toilets is a circle. No other 
things beside these are written. Obviously these two geometrical forms 
suffice and are fully understood by the whole population. 

In general dimension we talk about spheres. A sphere in d dimen- 
sions has its volume, V, proportional to r¢ and its surface, A, to r?1, 
so that V x AY(4-)), 

The area of sphere in three dimensions is A = 4ar* and the vol- 
ume V = 4nr3.. When V = 1, r = —(-3/m)'? (22/3, (3/m)' /22/3, 
or (-1)2/3 (3/m)!/? /22/3 which are numerically —0.310175 — 0.537239, 
0.62035, —0.310175, —0.310175 + 0.5372392 in that order. Therefore a = 
Alyai = 67/3,\/r = 4.83598 


2 


t Written while the author was at Chemical Engineering, UMIST, Su- 
pervisor, Professor Graham Arthur Davies, The University of Manch- 
ester, 7*” October 2002 
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When V = 8,r = (6/7)'/? and therefore A = 462/37!/3 = 19.3439. In 
order to find a, the surface area per unit volume, one divide the area 
by V?/5, in other words a = V'/9- A/V = A/V?/*. 

The same is true for other polyhedra. For example in a tetrahedron 
where z is the length of the side, the vertices can be 


gr V3 x V3 1 [93 
(0,0, 0), (x, 0,0), (5 F.0) ? (5 ene 2 =] 


When V = 1, one can obtain x by solving the equation 
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This gives x = 2 (38) = 2.0409 as the only real positive answer. When 
az is doubled, V increases from 1 to 8, which means that one would be 
dividing A by V3 to obtain a. 

The perimeter of a triangle is s = 3d. and the area A = (V3/4)I’. 
When A = 1, d = +2/3!/4 = £1.51967. Therefore s = 4.55901. 

Vertices shared by two cells make up a common face between them. 
Two way have been tried for finding the edges. The first one was by 
looking at all neighbouring cells of every cell in turn three at a time. 
The edges are then made up of those vertices that are common among 
these three cells. Only those edges which have exactly two vertices are 
considered. They are called good edges as contrasted with edges on the 
boundary. This is a much longer way than the second one, which is to 
consider vertices common to any two faces of a cell. Similar to the first 
case, such vertices forms a good edge if and only if there are only two 
of them. The two methods above give exactly the same list of edges, 
so they confirm each other. It has been tested that all edges having 
more than two vertices are boundary ones, that is they have at least 
one vertex outside the boundary of the unit cube considered. 











By drawing some of the cells as a solid using fill command it has 
been tested that the result from convhull covers the entire cell surface. 
This confirms the step where areas are calculated. 

The hexagon or honeycomb is perhaps the pattern which is most 
frequently found in nature. Even though the world we live in is three- 
dimensional, cells normally divide and spread in two dimensions in the 
form of layers. Moreover, they are packed in these layers in patterns 
which most often resemble the honeycomb (cf Williams and Bjerknes, 
1972). 
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An octagon is an eight-sided polygon. It is the shape of the cross 
section of every chimney in the mills built in Manchester during its 
industrial era of the nineteenth century, as well as that of the turrets 
in the Main Building of UMIST. Perhaps one of the reasons for its pop- 
ularity is that it looks strong while having the style of a good taste. 
May be the reason why it looks strong is that it possesses eight axes of 
symmetry, on top of another symmetry around the origin. 

There are nine regular polyhedra. Among these are five regular 
convex solids known to the ancient Greek called Platonic polyhedra. 
They are tetrahedron, cube, dodecahedron, octahedron, and icosahe- 
dron. They have regular congruent faces and regular polyhedral angle 
vertices. Their face angles and their dihedral angles at every vertex 
are equal. The other four regular polyhedra have only been discovered 
much later and are not convex. They are called the Kepler-Poinsot poly- 
hedra and are nonconvex. The small stellated dodecahedron and the 
great stellated dodecahedron were found by Kepler (1571-1630). The 
great icosahedron and the great dodecahedron were found by Poinsot 
(1777-1859). The small stellated dodecahedron and the great dodeca- 
hedron do not satisfy Euler’s equation. The process of creating it by 
extending nonadjacent faces until they meet is called stellating. There 
are also polyhedra called quasi-regular. 

The semi-regular polyhedra are called the Archimedean polyhe- 
dra. Here all faces are regular polygons but not all are of the same 
kind. Every vertex is congruent to all others. They comprise of an 
infinite group of prisms, an infinite group of antiprisms or prismoid, 
and another thirteen polyhedra. Each prism or prismoid is made up of 
two regular polygons on parallel planes where the vertices are aligned 
in the former case or shifted half way to the next neighbouring vertices 
in the latter case. Each vertex in prisms is joined to a corresponding 
vertex of the opposite polygon, while in prismoid it is joined to two 
corresponding vertices. All faces of an Archimedean solid are regular 
and all its polyhedral angle vertices congruent. 

On the other hand the Archimedean duals have the property that 
all their faces are congruent to one another and all their polyhedral 
angles regular. These solids are important in crystallography. They are 
vertically regular and include an infinite group of dipyramids, an infi- 
nite group of trapezohedra, and additionally thirteen other polyhedra. 

The surface area per unit volume a of a solid can be computed from 
the actual volume V and the actual surface area A as a = V‘/3A/V = 
AV~?/8, 
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Polygon Ne A s (numerical)| 
Triangle 3 (v3/4)d 4.55901 
Square 4 @& 4 
Pentagon 5 5(1+¥5)d2/ [420 —2/5)'/2] 3.8119 
Hexagon 6 33d? /2 3.7224 
Heptagon 7 = (7/4)d* tan(57/14) 3.6721 
Octagon 8  2d* tan(37/8) 3.6407 
Nonagon 9 (9/4)d? tan(77/18) 3.6198 
1/2 
Decagon 10 5 [6 +5) /2] @/(-1+ V5) 0.3605 
Undecagon 11 (11/4)d* tan(97/22) 3.5944 
Dodecagon 12 3(2+¥V3)d? 3.5863 


Table 1 Perimeter per unit area of n-gons. 


The tetrahedron is self-dual. The octahedron is dual to the cube 
while the dodecahedron the icosahedron. 

The nearest neighbour and minimum spanning tree have been ap- 
plied to the problem of taxonomy in botany. Clayton (1972), work- 
ing on the characters of plants to manually classify them (eg Clayton, 
1970) with the use of only the binary dendrogram and trial and error, 
adopted a numerical method which finds the minimum spanning tree 
in a multi-dimensional character space. Since taxonomy can be consid- 
ered as a kind of dictionary, it is possible to apply a similar approach 
to machine translation and the compilation of dictionaries. 

The icosahedron is common shape found among viruses. 

The polyhedra from Figure 1 to 3 are semi-regular. 





(a) (b) 
Figure 1 (a) Truncated tetrahedron, triakistetrahedron, 2 3|3. (b) Octa- 
hemioctahedron, octahemioctacron, 3/2 3|3 
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Figure 2 (a) Tetrahemihexahedron, tetrahemihexacron, 3/2 3|2. (b) Trun- 
cated octahedron, tetrakishexahedron, 2 4|3 


(a) (b) 





(a) (b) 


Figure 3 (a) Truncated cube, triakisoctrahedron, 2 3|4. (b) Rhombicubocta- 
hedron, deltoidal icositetrahedron, 3 4|2 


Polyhedra in Figure 4 are snub polyhedra. 
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(b) 


Figure 4 (a) Pentagrammic crossed antiprism, pentagrammic concave del- 
tohedron, |225/3. (b) Pentagrammic antiprism, pentagrammic deltohedron, 
|225/2 


The surface of Fullerine is made up of pentagons six-sided fig- 
ures. Its shape represents that of the geodesic domes developed by 
Buckminster Fuller, and hence the name Fullerine. The latter may ei- 
ther be hexagons or figures all the six sides in each one of which form 
two sets of three sides having an equal length. The simplest Fuller- 
ine, the carbon-60 molecule, has the same shape as that of a football 
and a handball. With some thought the reason for this is not difficult 
to see. With its thirty-two faces it closely resemble the sphere. Also 
the two different shapes of all its components are symmetrically dis- 
tributed and therefore enable colouring with only two different colours, 
namely one for each of the two shapes. To see how this helps, suppose 
one made a football in the shape of a bloated dodecahedron. Then it 
would be impossible to colour it using more than one colour at the 
same time of giving it a symmetrical appearance when viewed from 
more than a few directions. With the Fullerine shape and the colour- 
ing scheme mentioned, however, the football looks symmetrical when 
viewed from 54 different directions symmetrically distributed around 
it. These directions corresponds to those when one looks at it in the 
direction perpendicular to the centre of each of its faces and when in 
the direction through the middle of each of the 22 edges lying between 
two hexagonal faces. 

Making polyhedron models is an educating experience. Contrary 
to the general believe that you need to make an accurate drawing for 
the required parts (Wenninger, 1971), this needs not be so. Examples 
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of this are the origami models of polyhedra where complex polyhedron 
structures are made from interlocking pieces each of which is made by 
folding a piece of paper of a rectangular or square shape. 

A set of elements with the sum and the product of any two el- 
ements defined is a commutative ring if under these two operations 
it satisfies the following postulates: closure, uniqueness, commutative, 
associative, and distributive laws, identity (zero and unity), and addi- 
tive inverse. An integral domain is an ordered domain if its positive 
elements satisfy the laws of addition, multiplication, and trichotomy. 
A subset of an ordered domain is well-ordered if every nonempty 
subset of it contains a smallest member. a|b means that .b is divisi- 
ble by a. The Euclidean algorithm or division algorithm states that 
a = bg+7r,0 <r < b. Two integers are relatively prime if their only 
common divisors are +1. a = b(mod™m) if and only if m|(a — 6b). The 
commutative ring Z2 is the properties of multiplication and addition of 
even (0) and odd (1) numbers. 





+ 01 01 

0 01 0 00 

1 1 0 101 

The following is Zs. 

+0123 4 - 0123 4 
00123 4 000 00 0 
1123 4 0 1012 3 4 
223 4 0 1 202 4 1 3 
3 3 40 1 2 3 03 1 4 2 
4 40 12 3 4043 2 1 


There is a close link between geometry and algebra. Geometri- 
cal surfaces can be described as algebraical equations. For example, 
for circles and polygons the equations are binary quadratic, while for 
spheres and polyhedra they are ternary quadratic. Even one-sided sur- 
faces can be described algebraically. The equation of Klein bottle, when 
deformed into a sphere with two circles removed and replaced by two 
cross-caps, is a quartic equation 


a? (a ze y*) (b* _ x = y’) = 2*(a*a* aT by’), 
while the Steiner surface is also a quartic one 
ye + 2a? + 07y? + 2yz = 0. 
Two surfaces is homomorphic to each other if it is possible to con- 
tinuously transform one into the other. All convex polyhedra are homo- 
morphic to a sphere. The Steiner surface is homomorphic to the hepta- 


hedron, which is an Archimedean polyhedron with diametral plane. 
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In the plane, a second-degree equation gives either two straight 
lines, a circle, an ellipse, a parabola, or a hyperbola. In space, it can 
give two planes, cylinders and cones (circular, elliptic, parabolic, or 
hyperbolic), sphere, spheroid, ellipsoid, two hyperboloids, and (elliptic 
or hyperbolic) paraboloid. 

Partition, tessellation and division of space are the same thing. In the 
context of set theory, a partition of set X is a family of sets Aj, Az, ..., 
A, which are subsets of X, such that A; £4 0; A; A; = 0; U; Ai = X, 
where i, j = 1, 2,...,k andi# j. (cf Berge, 1958) A further condition 
that makes any tessellation a Voronoi one is that, for all 7 there exists 
a unique point a; within A; such that every point in A; is closer to a; 
than to any other a;, j # i. 

Voronoi tessellation in three dimensions can be constructed by 
imagining each region as a spherical cell growing outwards to meet 
neighbouring cells and continue growing to fill the gaps. The centre 
of each sphere is a unique nucleus point of the region such that it is 
closest to any point belonging to that region than any nuclei points. If 
the rate of growth is the same from every cell, the resulting partitions 
will be planes which can be described by ternary quadratic equations. 
However, if this rate differs from one cell to another, the partitions will 
be curved surfaces and the result is a non-Voronoi tessellation. It is 
possible to impose a constraint of minimum distance between neigh- 
bouring nuclei. Such cases can be looked at as spheres of an equal 
nonzero radius expanding away from nucleus centre points. If the radii 
differ from one sphere to another, or if some nonspherical solids are 
used instead of spheres, the tessellation obtained will be non-Voronoi. 


Consider the case where all spheres are of equal size. If these 
spheres already touch their neighbours before the expanding starts, the 
case is that of packed spheres expanded to form a Voronoi tessella- 
tion. There are two types of close-packing: cubic (face-centred) and 
hexagonal. In both cases each sphere has twelve neighbours. Both 
cases have the same density, which is =%5. The Voronoi regions pro- 
duced from the cubic case are rhombic dodecahedra and the faces are 
rhombuses. In the case of hexagonal close-packing, the correspond- 
ing regions are trapezo-rhombic dodecahedra and the faces are either 
rhombics or trapezia. Where the spheres meet with their three neigh- 
bours in the layer above and their three neighbours in the layer below, 
the faces are rhombics. Where they meet with the six neighbours on 
the same layer they are trapezia. 

For geometrical calculation, an example of a definitive book is that 
written by Salmon (1912). 


The gamma function, 
r= fo ete tat, Re(2) >0, (1) 
0 


12 January 2005 Vaen Sryayudhya, Editor 


Vol. 2, No. 1 Tyabandha Journal of Arts and Science 


got its name from Legendre and is known as the Euler gamma function 
or simply the second Euler function. The formula I'(z +1) = zI'(z) = 2! 
recursively calculates the gamma function from, for instance, T'(1/5) = 
4.5908, T(1/4) = 3.6256, [(1/3) = 2.6789, T(2/5) & 2.2182, T(1/2) = 
fa & 1.7725, T(3/5) & 1.4892, T(2/3) & 1.3541, T(3/4) & 1.2254, and 
T(4/5) = 1.1642. The Stirling’s formula was found by de Moivre which 
approximates the gamma function. The gamma function expansions is 


T(@ +1) = lim yeast (2) 


k00 (@ +1)(@ +2)---(@ +k)’ 
and the gamma function of negative numbers can be obtained from 


TT 





I(-z) = —————_-. 3 
ao zI(z) sinnz (3) 
The incomplete gamma function is 
T(z,2) = i evi / edi: (4) 
0 x 
and the normalised or regularised incomplete gamma function is 
T(z,2) 
T(z) 
Statistics 


Poisson distribution, defined by p(x,A) = [A*e~*/z!] Ijo,00), is the 
binomial distribution, p(x,n, p) = "C,6"(1 — @)"~*lio. nj, when n goes to 
infinity, 6 goes to zero, while n@ = ». Here @ is the probability of 
success of each trial. It is used when counting the number of occur- 
rences of a random event. Analogously Poisson point process, which 
has p(x = n(v)) = [Ajvle7*!”'* /z!] Ifo,00),, is the binomial point process, 
p(z =n(v)) = "C6" (1 — 8)" Tio nj, when the volume V goes to infinity, 
while n/|V| = ». Here 6 = |v|/|V| is the probability of points within V 
being placed in v C V c R¢, and \ the density or intensity of points. 
Therefore the density of point of a Poisson point process is constant by 
definition. A point process is a procedure which generates points on a 
domain within a space of d dimensions. 

The Poisson point process thus derived has the properties that 
0 < Prw)y=o0 < 1 for 0 < |u| < o, limjjop(n(v) > 1) = 0, n(v,) mu- 
tually independent and n(U,, v1) = >, n(vi) when A; are disjoint, and 
lim|e|-s0 [p(n(v) > 1)/p(n(v) = 1)] = 1. 
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The weighted mean of a group of data is x = )0; fiai/n and the 
weighted variance is o* = )), fi(aj — x)?/n, where f; is the occurrence 
frequency of x; and >, f; =n. Likewise the r'"-moment around the av- 
erage is m, = 0; fi(e—x)"/n, while the r'"-moment around the origin is 
m, = >>; fiz" /n (cf Spiegel, 1975). Some relations among these various 
moments are m, = 0, m2 = m4 — m,, mg = ms — 3mm, + 2m!?, and 
m4 = mi, — 4mm + 6m! mi, — 3m". 

The variance when normalised by n — 1 gives the best unbiased 
estimated variance if the sample has a normal distribution. On the 
other hand the variance which is normalised by n is identical with the 
second moment of the sample about its mean. 


Phase transition 


In the Ising model each spin has two possible states, that is up and 
down, and the Hamiltonian is H = Jo >); ; ci0; where the summa- 
tion is over the nearest neighbours. Since it has been exactly solved, 
the Ising model provides a good model for the understanding of phase 
transition. This model can represent the transition from ferro- to param- 
agnetic at the critical temperature where the correlation length becomes 
infinite. Characteristic to the Ising model is the peak in the specific heat 
at the critical temperature. 

The two-dimensional xy model is a model of spins confined to a 
plane, the Hamiltonian of which is H = Jo >7<; ; cos(@; — 9;). This 
model can represent the superconducting and the superfluid films. For 
this model there is no phase transition showing long-range ordering. 
One example is the two-dimensional Coulomb gas model where the 
vortex-antivortex pairs, which are bound to each other at low temper- 
ature, increases in number as the temperature increases and become 
separated at the KT temperature that marks the phase transition. 

It had been generally believed that no phase transition can ex- 
ist for the xy model when Kosterlitz et al (1973) showed that there is 
another kind of phase transition, arisen from the topological excita- 
tion of vortex-antivortex pairs instead of from the long-range order- 
ing found in a spontaneous magnetisation. They consider the two- 
dimensional model of gas with charges tq where the interaction po- 
tential is U(|r; — 1;|) is —2qq; In |(ti — 1j)/to| + 2u when r > ro, and 0 
when r < ro. The problem is reduced to that of solving an equation 
of the form (dy/dz) = —e~*¥. The application mentioned there is in 
the xy model of magnetism, the solid-liquid transition, and the neutral 
superfluid, but not in a superconductor and a Heisenberg ferromagnet. 

The frustrated xy model, the Hamiltonian of which is 


H = Jy S~ cos(6; — 6; — Aaj) 


<i,j> 
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occurs when a magnetic field is applied perpendicular to the two- 
dimensional plane of the xy model. The frustration parameter, f = 
@/Phio, is a measure of the average external magnetic flux. When 
f = 1/2 the model is called the fully frustrated xy model. The local 
chirality, m(r;) = x >5 (0; — 6; — Aij), which describes the property of 
the ground state, where it can either be +1/ovr2 or —1/ovr2. The net- 
work configuration at T < T, is that of a draught board, and has Z» 
symmetry. This regularity is broken by the formation of domain walls 
in an Ising phase transition at T>. 


Random processes 


A synonym to random is stochastic (cf Miles, 1972). Any algorithm 
which employs a random element is called Monte Carlo. Random pro- 
cesses can have various types of distribution. The beta distribution has 
a probability density function fx, ,(«) = oe 0 <a <1, where 
a >0Oand 6 >0 are shape parameters, and B(a, 8) is the beta function. 
There are three types of shape; the bridge shape has a > 1 and f > 1, 
the J shape a <1land 6>1, ora>1and 6 <1, and the U shape a < 1 
and 8 <1. 

The 200 generators used in Figure 5 are randomly chosen with beta 
distribution with the shape parameters a = 2.7 and £ = 3, that is bridge 
shape. 




















Both z and y in Figure 6 have J-shaped distribution with the shape 
parameters a = 2 and 6 = 0.8. The density is unbounded at x = 1 and 
at y = 1 because 6 < 0 for both. 
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The shape parameters in Figure 7 are a = 0.5 and 6 = 0.3, that is 
U shape. The density is unbounded at x = 0, 1 and at y = 0, 1 because 
a is also less than zero.: 




















Structures in nature 


Prusinkiewicz and Lindenmayer (1990) model the structures in 
plants after having briefly discussed about the difference between the 
Chomsky grammars and the L-system, both of which being a mean for 
doing string rewriting, but the latter is based on the Turtle geometry 
which makes it convenient for the geometric rewriting of fractals. The 
states of a turtle consists of its position coordinates and the direction in 
which it is facing. Meinhardt (1995) models the patterns on sea shells 
by using mathematical based on the partial differential equations gov- 
erning the system of activator, inhibitor, and substrate. Starting from 
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a homogeneous initial condition, small deviations therein undergo a 
positive feedback and therefore increase. Activator catalyses both the 
production of itself and that of its inhibitor. The latter acts as a negative 
feedback which limits and makes the reaction local. 


Random tissue in three dimensions has four edges, six faces and 
four cells meeting at each vertex. It is thus surrounded by four cell 
nuclei, as well as by six bonds forming a tetrahedral cage. The four 
edges meeting at a vertex resemble a caltrop, and similarly the lines 
from it to the four nuclei. If a vertex had more than four edges it would 
have been structurally unstable, because then it can be split into two 
normal vertices by an infinitesimal deformation. Continuous random 
networks, for example models of covalent glasses like vitreous silica, 
are excluded from this restricted class since theirs may be more than 
three non-planar faces to an edge even if they still have four edges to a 
vertex (Revier, 1982). 

Unlike crystallography, the ideal random structure is by no means 
unique because it is the solution of a statistical problem. There are, 
however, certain geometrical and topological invariances, the most fa- 
mous of which is possi bly the Euler’s theorem. In two dimensions this 
theorem states that f—e+v = x, where x is an Euler-Poincaré character- 
istic and integer of order one, x being for instance 1 and 2 respectively 
for plane and sphere; in three dimensions it is f-—e+v = 2. The valence 
relations, So nf, = 2e = 3v, hold for 2-d and 3-d alike. 

Stumbling upon some observations, Theorem’s 4 and 5 resulted. 
These two theorems help explain together with Algorithm & on page 
X, the valence relations. Theorem 1 is also another product obtainable 
from applying both the Euler’s theorem and the valence relations (cf 
Prause, 2000). 


Theorem 1. The average number of edges per polygon in a large pat- 
tern is six.: 


Proof.: From Euler’s theorem, f — e+ v= 1, and the valence relations, 
Yndn = 2e = 3v, it follows by applying the latter to the former that 
f —dnf/2+ donf/3 =1. Since i is the average number of edges per 
face, it follows that f—nf/2+7f/3 = 1. Then f(1—n/6) = 1, and 
consequently 1—1/f = 7/6. As the network becomes large, f becomes 
infinite and as a result 7 = 6. o 

It is interesting to note that Theorem 1 posts no restrictions on 
whether the network is random or regular. The only assumption made 
is that three and only three edges meet at each vertex. One is almost 
tempted to say that, as the size of a network becomes infinite, nature 
somehow follows this theorem and make sure that each of the polygons 
has six edges on average. 


When I translated the three papers by Voronoi (Tiyapan, 2001) I 
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used the term vertex to mean a vertex of a specific polygon or polyhe- 
dron. For any vertex, I used the term vertice, and for more than one 
vertex vertices. It turns out that I am not the only one who concerns 
himself with the word. Moore and Angell (1993), for instance, use apex 
for a single corner, vertex to mean any point in a tessellation in two or 
three dimensions where its apices meet. 

Rivulets flowing inside a pack bed have been studied by several 
authors (cf Porter, 1968). They are said to flow independently of each 
other, with no mixing among one another. Diffusion theory has been 
used to treat a random walk process. 

The van der Waals force in combination with double layer repul- 
sion play an important part in the study of filter and particle move- 
ments in porous media. Both are electrical forces and can be used to 
explain the particle capturing mechanism. 

Percolation is related to chemical engineering contexts (cf Mohanty 
etal, 1982). But the necessary background in stochastic processes for the 
purpose of simulation has already been described earlier, for example 
the modelling of colmatage, the retaining of particles suspended in a 
fluid flowing through a porous medium (Litwiniszyn, 1963 and 1967). 
In general, however, few authors in all engineering fields relate their 
works directly to the percolation theory. The majority of studies in this 
area are based on dynamics and fluid dynamics theories (cf Mulder and 
Gimbel, 1990; Kock and Judd, 1965). 

The reason for the lack of percolation material in engineering lit- 
erature may be because the percolation process generally works behind 
the scene and only shows itself as critical phenomena. Most engineer- 
ing studies are concerned with things under some operational condi- 
tion, within the range of which percolation seems to be absent. By 
contrast, in Physics where extreme conditions are considered, there is 
an enormous and increasing amount of publications which are directly 
under the topic of percolation. But this by no means means that within 
more pacific ranges no places exist for percolation. In fact it is precisely 
this lack of mentioning in the literature that induces one to this kind 
of research. Because our idea of criticality is by tradition closely linked 
to the idea of time, percolation seems to be present only when there 
are instantaneous changes. But if in a steep s-curve we only rotate the 
axis clockwise by 7/2 radian, such that to make the time axis vertical 
instead of horizontal, then we will see that in place of one critical point 
in the middle of the graph connecting two different levels, there are 
now two critical phenomena on both sides, one on each side, and in the 
middle a flat region where time hardly changes. Looking at it this way 
percolation seems to be a symmetry between time and space. In physi- 
cal systems spaces percolates, but in the dual world where criticality is 
continuity it is the time instead which percolates. 
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Having said that, without rotating the time-space axes backwards 
and forwards too often one should still find ways to investigate what a 
pacified percolation does. In this regard, the study of economics seems 
to be an ideal place to start, if simply because one knows there exists 
such thing as hyperinflation but one never wants to study that when 
it happens. This automatically forces the researcher to find ways of 
doing researches which would not ruin his pocket or put his life in 
jeopardy. Another ground with a good prospect is in traffic congestion, 
even though its worst effect is not yet devastating, apart from what it 
sometimes does to the economy. 

With these digressions in mind, if we now turn our thought at this 
point to our chemical engineering studies, we will not fail to see how 
ideal filtration fits the no-ruins requirement. For here we have a process 
which percolates routinely, often needs to be backflushed, but where 
the effect it produce is probably that of giving engineers a headache 
and, at its worst, putting a decent company out of business. But even 
with this advantage, I still believe that the study of filtration should not 
concentrate only on the fouling of filters, but should try to understand 
both the percolated and nonpercolated situations, preferably the latter 
for the lack of it in literature, and to connect what happens in a working 
filter to what happens, or does not happen, in a fouled one. 

The revolutionising discovery made by F. August Kekulé (Kekulé, 
1865, cited in Wotiz, 1993) that benzene has cyclic nature gave rise to 
the structural theory of organic chemistry. 


Voronoi tessellation 


In a puzzle of Figure 8 there are five rooms with doors in the 
position as shown. The problem is whether one can walk through 
every door only once and the answer according to the graph theory is 
no, because there are more than two rooms which has an odd number of 
doors. The proof of Theorem 2 was from Komsan Bajaravanijy around 
1989. From this, when one plays a puzzle like that of Figure 8 one 
always starts off from a room with an odd number of doors and ends 
in another such room. This is the same thing as saying that one starts 
and ends outside rooms with an even number of doors. Therefore the 
number of the latter is of no consequence, but that of the former is 
crucial for the existence of a solution and must never be any number 
other than two. 

In Figure 8 there are three rooms with five doors, two with four, 
and one with nine. There are here four rooms with an odd number of 
doors. Starting off from one of these four one can only end up in one 
of the other three, which leaves the remaining two rooms unaccounted 
for. In other words at least two doors will necessarily remain unvisited. 
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Theorem 2. A travel along a network can only starts and ends at 
nodes which have an odd coordination number. 


The proof of Theorem 2. 


Proof.: Looking at Figure 9 if one starts from inside a room with an 
odd number of doors, one always ends up outside it. On the other 
hand one always ends up inside a room with an even number of doors 
if one starts of from it. In the second picture such a room is all the area 
outside the circle. 
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The following Corollaries 2[1], 2[2] and 2[3] assume nondegeneracy 
of the Voronoi network. Such a path as mentioned in these corollaries 
is also called self-avoiding. 


Corollary 2[1]. There can be no path which traverses all edges of a 
Voronoi graph only once. 


Proof.: This follows from Theorem 2 because a two dimensional Voro- 
noi network has a coordination number three. o 


Corollary 2[2]._ Ona three dimensional Voronoi structure there always 
exists a path that runs through every edge once and only once. 


Proof.: This also readily follows from Theorem 2 since a three dimen- 
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sional Voronoi network has a coordination number of four. a 


Corollary 2[3]. Take any Voronoi cell of the three dimensional net- 
work, it is impossible to walk through all its edges without repeating 
some of them. 


Proof.: Again from Theorem 2 and because the surface of a Voronoi 
polyhedron is a two dimensional network of polygons which has the 
coordination number three. o 

Jerauld et al (1984) compared the Voronoi, with the triangular net- 
works and found that the bond percolation probability of the former is 
4.3% or 0.015 less than that of the latter, small site clusters more, and 
small bond clusters less likely. 

When a square lattice was fed to voronoi and delaunay in Matlab, by 
the programme degen.m, there were error messages saying that points 
were collinear and possibly triangulation is incorrect. This case, Figure 
10, is degenerative. 


Voronoi and Delaunay, degenerate case 
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Distorted honeycomb, one third shift Distorted honeycomb, two third shift 
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Figure 11 The honeycomb or hexagonal lattice whose alternate y-plane has 
been shifted (a) one-third, and (b) two-third respectively. Triangulation is 
shown with thinner lines. The programme used is honey .m 
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(a) a © 


Figure 12 (a) The hexagonal lattice, (b) its covering lattice (Kagome), and (c) 
the covering lattice of its covering lattice (i.e. covering of Kagome). (cover .m, 
covers.m,and coverss.m) 


If we indicate by C(x) the n*-order covering lattice of a lat- 
tice x, then the first picture is Hexagonal, the second one Kagome 
or C}(Hexagonal) and the third one C?(Hexagonal) or in other words 
Cl(Kagome). Figure 13 is the next iteration, a C3(Hexagonal) or 
C?(Kagome). 
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Covering of covering of Kagome 




















Now let us look at the Voronoi graph and its covering lattices. 
Pictures in Figure 14 are drawn by first creating and cropping a Voronoi 
graph with the help of the programme crop.m, then use the recursive 
procedure described above to find up to the third covering lattice. 


Voronoi graph 
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Second covering of Voronoi graph Third covering of Voronoi graph 
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Figure 14 (a) Voronoi graph (V.g.), (b) Cl(V.g.), (c) C2(V.g.), 
(d) Cy (V.g-). 


Notice that the covering lattices retain the skeleton structure of 
the original Voronoi graph. Those of higher orders represent closer the 
structures of nature where walls have thickness. 


Quadratic equations 


Quadratic equations are equations of binary quadratic forms. 
Around 400 BC Barbilonia had algorithmic equivalences of quadratic 
equations which are based on the method of completing the square 
and where all answers are unsigned, i.e. positive, lengths. Because 
there was no notion for zero, Diophantus considered three types of 
quadratic equations axz* + bx = c, ax* = be +c, and az? +. ¢ = ba. Eu- 
clid, circa 300 BC, used geometric equivalences or quadratic equations 
whose roots are also lengths. Brahmagupta allowed negative quantities, 
which he called debts, and used abbreviations for the unknown. Al- 
Khwarizmi classified quadratics into six types, namely squares equals 
roots, squares equals numbers, roots equal numbers, squares and roots 
equal number, squares and numbers equal roots, and roots and num- 
bers equal squares. In his book Liber embadorum, published in 1145, 
Abraham bar Hiyya Ha-Nasi (aka Savasorda) gives the complete solu- 
tion of quadratic equations. Luca Pacioli published Summa de arithmetica, 
geometrica, proportioni et proportionalita (or Summa) in 1494. He also ap- 
plied quadratic methods to quartics of the form x* = a + bx”. Scipione 
del Ferro solved the cubic equations of the form 2° + mz = n. 


Quadratic forms 
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The theory of quadratic forms and the theory of matrix are insep- 
arable though the history of these two subjects are somewhat fragmen- 
tary. A bilinear form in the sets x; and yj, i = 1...n, is pare L4Yj, OL 
xT Ay. If x; = y; for all i, then the form is quadratic in x;. In other 


words, a quadratic form is a general expression which contains second 
order terms. 


Voronoi statistics 


In two dimensions the statistical descriptions are as follows. 

The curve is y = 0.62/n°’.: First the edge lengths are normalised 
by the edge length of the equivalent or characteristic square . A equiv- 
alent square is defined as the square figure whose area is equal to the 
area of the polygon in question, here a Voronoi polygon. Then find the 
average of the edge lengths in each cell. 





Minimum normalised edge length 
roy ry 














f f f 
10 10' 10° 10° 10° 10° 
Number of cells 


The curve is y = nis)" + 0.7.: Of these cell-averaged nor- 
malised edge length obtained from simulations on various sizes of net- 
works the minimum values are plotted in Figure 15, the maximum in 
Figure 16, and the expected value in Figure 8. Note that the last quan- 


tity is the average over the whole structure of all the averages obtained 
one from each cell. 
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Maximum expected normalised edge length 
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The curve is y = 10°?/!98* — 0.4.: To summarise, as the networks 
gets larger its minimum, maximum, and mean of the edge lengths when 
compared with the characteristic length approach constant values. The 
characteristic length is defined as 1 = [}7j", A; /nj'/ * where A; is the 
area of the i'* polygon and n is the number of cells. 





Mean expected normalised edge length 
roy 
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That the three values mentioned become constant may not seem 
obvious by the look of Figures 15 to 17 because the scale used there 
is a logarithmo-logarithmic scale, not a Euclidean one. These figures 
emphasise the smaller ranges of size. Figure 18 below on the other hand 
is plotted using a normal scale which enables one to see the asymptotic 
effect more clearly. Here the figures (a), (b), and (c) are respectively 
Figure 15, 16, and manel. Let the term representative stands for ‘of the 
cell-average normalised’, and length means ‘edge length’ in this context. 
Then the minimum representative length approaches 0.15 from Figure 18 (a), 
the maximum representative length is ever increasing, seemingly by a power 
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law of approximately 0.5, while mean representative length approaches the 
value of 0.65. Properties of the Voronoi tessellation can be divided 
into individual and collective properties. With this in mind the term 
length above represents a property, representative means individual, and 
minimum, maximum and mean show the collective attributes. 
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This o?(N¢(E(I-))) increases very slowly with the increasing sizes 


of the networks. The curve shown has the equation 


28 


y= [2/8 x 108)]"| +0.05 
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The number of vertices per cell increases dramatically as one goes 
up the dimension ladder. The programme used contains the essential 
part of the code which produces Figure 20. The straight line shown 
is 0.093(4 + e)”. Notice the trend towards a greater rate of increase at 
dimensions higher than the maximum six shown. Only Voronoi cells 
which lies within the original domain and the vertices of these are 
considered. The same programme also gives Figure 21. 
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The expanse and the hypervolume of the Voronoi tessellation in- 
creases at an enormous rate, which results in the number of cells totally 
bounded within the original domain decreasing rapidly in Figure 21 as 
the dimension goes up. The effect at close to the zero ratio is empha- 
sised by using the log scale for the y-axis. Also with the logarithmic 
scale the polynomial estimation curves that give negative values of the 
ratio is automatically excluded. 
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Dimension 

One can fit the polynomial p(x) = piz” + pox" 1 +...+ pnt +pn4i to 
the data with a least square algorithm. If x is the vector containing the 
data, then the the (n + 1) coefficients of the estimated polynomial can 
be found from X = (x — E(x))/o(x). For data containing independent 
normal errors with a constant variance, the error bounds contain at 
least half of the predictions. The curve shown in Figure 21 is p(x) = 
—0.0482* + 0.06427 — 0.20327 — 0.459x + 0.263, the average value E(x) 
is 3.917, and the standard deviation o(x) is 1.412. The structure of 
the polynomial fit can be described using the Cholesky factor of the 
Vandermonde matrix 


—12.19 1.56 -6.08 -—-047 —3.17 
0 8.46 —0.45 447 —0A44 





R= 0 0 1.21 033 2.93 |, 
0 0 0 —1.63 0.30 
0 0 0 0 2.25 


the degree of freedom which is 19, and the norm of the residuals which 
is 0.034 in this case. 


Voronoi section 


Grouped area distribution 


Diagonal vertical section 
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min 8 6 2.2191 x 10-4 9.9722x 10-5 12 4 0.6030 
max 46 25 4.7015 x 1073 0.14364 69 5.52 2.5690 
lb 26.338 15.169 1.8975x10-% 1.8975 10-3 39.507 5.1712 1.5586 
o2 40.608 10.152 5.2144 10-7 4.2313 10-5 91.368 0.035983 0.1186 
o 6.3725 3.1862 7.2210x 10-4 6.5048 x 10-3 9.5587 0.18969 0.3444 
Mg 25.543 14.829 1.7572x 10-3 1.1575 x 107 38.315 5.1676 1.5169 
Mn 24.703 14.478 1.6026x 10-3 8.5548 x 10-4 37.054 5.1638 1.4700 
med 26 15 8125 x 10-3 1.1660 10-3 39 5.2 1.5706 
mad 5.0068 2.5034 5.6975x10-4 1.4250x 10-3 7.5103 0.14465 0.2715 
M? 40.531 10.133 5.2045x 10-7 4.2233x 10-5 91.195 0.035915 0.1184 
Me 72 9 2.5644x10-19 5,.4435x10-° 243 —8.3404 x 1073 -0.0060 
M* 5088.9 318.06 1.0467x10-'2 7.6653x 10-7 25763 7.9566x 10-3 0.041906 
K 3.0978 3.0978 3.8644 429.77 3.0978 6.1689 2.9916 





Table 2 Simulation uses rbox (1000 random points, seed 234985) and qhull (option v and 0); d = 3, n° = 527, 
n” = 6357, CPU time 6,466.99 sec for the counting of statistics, 270.56 sec for finding area of the faces, 4.47 sec for 


calculating cell volume and 2.8 sec for finding the number of edges. 
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Number of vertices and edges 


It has been observed from the simulations that in three dimensions 
cells always have vertices in even numbers and edges odd ones. This 
can be explained by the following theorems. 


Theorem 3. (cf Miles, 1972) In a simple three dimensional Voronoi tes- 
sellation, 3n% = 2né. 


Proof.: Pick any Voronoi cell of the tessellation. Suppose that it has 
n” vertices. Add up the number of edges connected to all vertices. Be- 
cause every cell is a simple polyhedron, there are exactly three edges 
connected to each vertex. The number of edges thus counted is there- 
fore 3n”. But each of the edges is connected to two vertices, so we have 
counted every one of them twice. Therefore, 


2n° = 3n”. 


This is the case for any cell, hence the theorem is proved. o 
This theorem gives rise the following two theorems. 


Theorem 4. The number of vertices of any cell within a simple three 
dimensional Voronoi tessellation is an even positive integer. 


Proof.: Observe that the term 2n¢ in the theorem above is divisible by 
2. This term is equal to 3nt, therefore the latter is also divisible by 2. 
Since 2 can not divide into 3, the only term left, n?, must be divisible 
by 2 and hence even number. o 


Theorem 5. The number of edges of any cell within a simple three 
dimensional Voronoi tessellation is a positive integer divisible by three. 


Proof.: With the same line of reasoning as above, observe that 3n® is 
divisible by 3. Therefore 2n€, and hence né§, is also divisible by 3. o 
Another proof for both the above theorems is the following. 
Proof.: For an equality to hold, both sides must have the same factors. 
By supposing an unknown common factor i and by cross-multiplying 
the coefficients on both sides, one obtains 2- (3-1) =3-(2-i), where i is 
a positive integer. Therefore nf = 3i and n? = 2%. In other words, né is 
divisible by three and n? is even. o 
The theorems above assume that every cells are simple. This can 
not be the case in real situation where edges have dimensions and 
rather represent tubes than one-dimensional lines. Such case is similar 
to the so-called degenerative case in a computational model of Voronoi 
tessellation where there exist vertices the number of edges connected to 
each of which exceeds four. Even in the degenerative case, one would 
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perhaps still expect a tendency for n® to be an even number and for n& 
to be divisible by three to hold. 

In nature there are things which have a tendency towards even 
numbers. The following graph shows the the abundance in cosmic ma- 
terials from the compilation by Cameron (1973). The year of publication 
of this paper is often misquoted as 1970, due to misprints in a footnote 
on the first page of the paper. Even the legendary Fred Hoyle has con- 
sistently practised this mistake and let it go uncorrected throughout his 
career. From what I have come across without an effort of searching, in 
two of his books and at least one of his papers, spanning the period of 
roughly forty years in total (cf Hoyle, 1977). Out of a sample of 278 of 
those papers which cite this work, this misprint resulted in 80% of the 
total number of errors in the year cited, which in turn amounts to 1.8% 
of the number of samples. 

Figure 23 shows The abundances of the chemical elements in the universe. 
They are assumed to be the same as those found in the primitive solar 
nebula, which have been deduced from data on abundances found in 
chondritic meteorites and those found in the Sun. All abundances are 
relative to that of Si which is taken to be 10°. Missing bars appear 
where the atomic numbers are unstable. Except for the atomic number 
1 of Hydrogen, which is the most universal element, all other elements 
with even atomic numbers are locally more abundant than those with 
near-by odd atomic numbers. 





Cosmic abundance 








0 10 20 30 40 50 60 70 80 90 100 
Atomic number 


Of interest are also the electrical resistivity and conductivity of 
solid matters. The conductivity o is by definition the reciprocal of the 
resistivity p. Some of the solids, particularly boron, carbon, silicon, sul- 
phur, germanium, selenium and tellurium, have a distinctively higher 
resistivity than the majority. Interestingly all of these, with only one 
exception of boron whose atomic number is five, are of an even atomic 
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number, which respectively from carbon are 6, 14, 16, 32, 34, and 52. 
This can be seen in Figure 24 the data of which are taken from Podesta 


(2002). 
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Figure 24 The electrical resistivity, (a), and conductivity, (b), of elements 
which are solid at the room temperature. 
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Ny (EC) 


Ne 
min 3 
max 11 
bb 5.8973 
o 1.8955 
a 1.3768 
Lg 5.742 
Ln 5.5904 
med 6 
mad 1.0736 
MM? 1.8912 
Me 1.6194 
M+ 12.438 


0.20716 
8.1072 

0.70626 
0.18607 
0.43136 
0.66036 
0.62941 
0.65709 
0.17995 
0.18565 
0.96371 
6.8376 

198.38 


NS (Gc) 


0.2749 
10.134 
1.0089 
0.2942 
0.54241 
0.94795 
0.89999 
0.96614 
0.24357 
0.29355 
1.7931 
15.785 
183.18 


Ny (Ac) 


0.058428 
22.307 
1 
1.4379 
1.1991 
0.79225 
0.61799 
0.84267 
0.49216 
1.4347 
22.703 
468.03 
227.39 


K 3.4775 


Table 3 Neighbour statistics. Here for the normalisation purpose, lf 


basis — 


0.046662, Qpasis = 0-18665, Apasis = 0-0021773. Simulation uses voronoin 
command in Matlab; d = 2, n° = 448, n” = 946, CPU time 1.19 seconds. 


Next simulation was done with d = 2, n° =3 to 49551. 

Figure 25 shows percentage of space covered by a Voronoi struc- 
ture. The number of cells is the total number of cells generated. The 
percent space covered is the volume of the structure after boundary 
cells, that is cells which extrude the unit volume boundary, have been 
excluded. The equation of the reference curve is y = —210/ log x + 120. 





Percent space covered 
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The reference line in Figure 26 is the linear equation y = 2x + 10. 
Boundary cells have been excluded. 
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The curve in Figure 27 is the result of curve fitting by cubic spline 
interpolation. 
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The characteristic length is the length of the side of the cubic struc- 
ture having the same number of cells and the same total volume as the 
Voronoi structure. The characteristic lengths in Figure 28 are shown as 
dots. The reference line is y = 0.7/a2°-®. 
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The characteristic area is the area of each square in the assembly 
the total volume and the number of cells of which are the same as 


those of the Voronoi graph. The reference curve shown in Figure 29 is 
y = 0.43/x°?. 
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Figure 30 Distribution of the number of edges per face. Each picture is an 
individual cell. The distribution shows the relative abundance or the number of 
faces (the vertical axes) having the number of edges as shown by the horizontal 
axes. The horizontal axis scales are positive integers starting from zero at the 
origin. Simulation uses rbox (500 random points, seed 893280) and ghull 
(option v and o); d = 3, n° = 231, n” = 3,107, CPU time 81.1 sec for finding 
face area, 2.22 sec for cell volume, 49.59 sec for counting edges and 0.03 sec for 
finding the number of edges and faces. 


The minimum number of edges for each face is three. This number 
is the same as the number of vertices of that face. From these figures 
most of the cells have at least one face with three edges. There are only 
19 cells (8.23 per cent) which does not have any three-edged face, and 
all of them have some four-edged faces. Therefore there is no cell with 
five as the minimum number of edges per face. The maximum number 
of edges per face is less clear-cut. There are twelve cells (5.19 per cent) 
with 11 as the maximum number of edges per face and two (0.87 per 
cent) with 12. 
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So0z Auunuv{ 
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N= 2 


2. 


UWOUTLOD UT X9}IAA dUO sea] ye SuTALY s][ao OM} Aue ‘spIOM IayjO UT *7U 
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Aopipy ‘vhypnhvhas uav, 





ne Re Afr ye a ne E(n§,,) o(n§..) 
min 10 7 1.0989 x 10-3 2.5790 x 10-4 0.0406 15 4.2857 0.6030 
max 44 24 1.0956 x 10-7 0.24840 4.2505 66 5.5000 2.7028 
U 25.974 14.987 4.3290 x 10-3 4.3290 x 107% 1.2515 38.961 5.1601 1.5450 
o> 40.852 10.213 3.3308 x 10-© 2.9419 x 10-4 0.3802 91.916 0.0372 0.1243 
o 6.3915 3.1958 1.8250 x 1073 1.7152 x 107? 0.6166 9.5873 0.1928 0.3526 
Ig 25.166 14.642 3.9633 x 10-3 1.8465 x 10-° 1.0816 37.749 5.1564 1.5030 
Mn 24.323 14.288 3.6026 x 10-3 1.2471 x 10-3 0.8241 36.484 5.1526 1.4580 
med 26 15 4.1467 x 10-3 1.5715 x 107% 1.1915 39 5.2000 1.5315 
mad 5.1532 2.5766 1.4064 x 1073 4.5864 x 10-3 0.4720 7.7298 0.1505 0.2794 
M2 40.675 10.169 3.3163 x 10-© 2.9291 x 10-4 0.3785 91.518 0.0370 0.1238 
M? 59.377 7.4222 5.7258 x 10-9 6.3817 x 10-5 0.2143 2.0040 x 10? -0.0071 0.0058 
M?! 4.5398 x 10° 2.8374 x 107 4.6684 x 10-1! 1.5395 x 10-5 0.7390 2.2983 x 104 0.0063 0.0463 
kK 2.7440 2.7440 4.2448 1.7943 x 102 5.1578 2.7440 4.5730 3.0230 


Figure 31 Statistics of a 231 cells Voronoi structure. 
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are neighbours to each other, and the number of neighbours around 
any cell in an infinite network is equal to the number of its faces. 

The line shown in Figure 32 is y = 0.2410e", where n is the di- 
mension of the network and is the horizontal axis; the coefficient of the 
exponential term is obtained by averaging over the averages in each 
dimension, each of which in turn comes from five batch runs. 


10° 





Number of vertices / number of nuclei 
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10° Ll L L Ll L L 
2 3 4 5 6 7 8 9 


Dimension 





The line shown in Figure 33, found manually by trial and error, 
has the equation y = 4.61 x 10~°(2 + e)”. In comparison, substituting 
the average cpu time for each dimension for y in the equation y = Ae” 
to obtain A and then find the average again over all dimensions results 
in A = 1.612 x 10-*. The programme which produces both Figure 33 
and 32 is given at the end. 
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The number of vertices of higher dimensions is investigated briefly 
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Dimension _ 4 





Ne 300 

Ny 6,577 
Ne 132 
min(n?) 26 
max(n”) 255 

ne 09.90 
(97, e 1,155.9 
on 33.999 
[g(n?) 04.41 
bn (n®) 98.389 
med(n®?) — 107.50 
mad(n?) 27.366 
m?(n2) 1,152.1 
m3(n?) 15,555 
m'(n?) 4.5760 x 10° 
K(n?) 3.4476 
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104 
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2,666.1 
2,571.4 
2,651.0 
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423.84 
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1,542.3 
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1.7365 x 10° 
4.7965 x 107 
8.5043 x 101° 
2.8202 


Table 5 Statistics of the number of vertices in 4, 5, 6, 8, 9, and 10 dimensions 
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Frequency of number of vertices per cell for a 4--d Voronoi of 300 cells, random seed 6566545, 
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Frequency of number of vertices per cell for a 5--d Voronoi of 300 cells, random seed 125980 
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Frequency of number of vertices per cell for a 6-—d Voronoi of 300 cells, random seed 467009 
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Frequency of number of vertices per cell for a 9--d Voronoi of 300 cells, random seed 75689 
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Figure 34 Distribution of vertices in (a) 4, (b) 5, (c) 6, and (d) 8, (e) 9, (f) 
10 dimensions 


To obtain the number of vertices per cell, n?, of a six-dimensional 
Voronoi structure of 1,000 cells I used a batch programme, for instance 
the one listed at the end. The Matlab macro that this programme refers 
to opens and reads from a file the number of vertices and then finds the 
statistical values, viz. min(n®) = 1,198; max(n?) = 9,923; #2 = 4,201.1; 
(o7,,)e = 1.4069 x 10°; o% = 1186.1; (iy)? = 4035.4; (rin)? = 3866.8; 
med(n?) = 4122.5; mad(n?) = 931.81; m?(n®) = 1.4055 x 10°; m3(n2) = 
1.0462 x 10°; m4(n®) = 7.7148 x 101%; and K(n®) = 3.9053. 

The following figure shows the frequency of the number of ver- 
tices. 


6--d Voronoi, 1,000 sites with random seed 3765098 
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There are three cells with 2729, 3213, 3246, 3421, 3490, 3606, 3938, 
4143, 4398, 4417, 4442 and 4970 vertices. There are two having 1854, 
2549, 2557, 2634, 2712, 2720, 2722, 2751, 2804, 2843, 2869, 2878, 2882, 
2920, 2931, 2957, 2973, 2996, 3013, 3090, 3260, 3322, 3329, 3343, 3366, 
3393, 3418, 3424, 3446, 3513, 3520, 3560, 3577, 3583, 3585, 3610, 3631, 
3677, 3714, 3732, 3753, 3767, 3770, 3819, 3835, 3861, 3907, 3919, 3924, 
3937, 3939, 4045, 4053, 4091, 4117, 4122, 4123, 4163, 4178, 4219, 4251, 
4261, 4269, 4272, 4298, 4317, 4327, 4355, 4461, 4470, 4473, 4474, 4542, 
4563, 4572, 4576, 4586, 4618, 4624, 4629, 4640, 4646, 4648, 4700, 4729, 
4792, 4810, 4851, 4942, 4975, 4984, 4985, 5058, 5064, 5097, 5161, 5195, 
5235, 5287, 5347, 5387, 5455, 5467, 5492, 5529, 5606, 5650, 5838, 5913, 
6112, 6193 and 6455 vertices. And there is only one cell with each of 
the following numbers of vertices, 1198, 1612, 1672, 1686, 1688, 1717, 
1727, 1787, 1827, 1864, 1906, 1915, 1921, 1947, 2016, 2052, 2056, 2060, 
2123, 2190, 2200, 2223, 2231, 2236, 2267, 2280, 2281, 2286, 2289, 2344, 
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2353, 2356, 2371, 2372, 2385, 2408, 2423, 2426, 2429, 2455, 2475, 2487, 
2499, 2500, 2503, 2504, 2508, 2513, 2542, 2547, 2551, 2560, 2561, 2565, 
2575, 2578, 2580, 2585, 2586, 2596, 2617, 2619, 2622, 2650, 2654, 2658, 
2664, 2669, 2673, 2686, 2692, 2694, 2704, 2708, 2716, 2718, 2723, 2732, 
2733, 2742, 2743, 2755, 2771, 2779, 2781, 2783, 2791, 2800, 2807, 2808, 
2811, 2820, 2835, 2837, 2853, 2858, 2864, 2867, 2870, 2873, 2875, 2880, 
2884, 2887, 2893, 2894, 2895, 2902, 2907, 2910, 2914, 2916, 2925, 2928, 
2934, 2949, 2955, 2956, 2960, 2967, 2974, 2978, 2983, 2985, 2989, 2993, 
3007, 3012, 3014, 3027, 3051, 3074, 3102, 3107, 3114, 3116, 3119, 3129, 
3130, 3136, 3137, 3141, 3148, 3152, 3164, 3173, 3176, 3180, 3182, 3186, 
3188, 3190, 3192, 3199, 3202, 3204, 3206, 3208, 3219, 3221, 3229, 3231, 
3237, 3239, 3244, 3245, 3255, 3256, 3259, 3267, 3270, 3271, 3274, 3275, 
3285, 3295, 3296, 3305, 3307, 3312, 3316, 3317, 3318, 3319, 3333, 3337, 
3344, 3351, 3352, 3354, 3357, 3360, 3370, 3373, 3374, 3378, 3384, 3390, 
3394, 3396, 3402, 3403, 3408, 3427, 3428, 3436, 3440, 3443, 3444, 3452, 
3453, 3454, 3486, 3487, 3492, 3499, 3502, 3505, 3506, 3509, 3519, 3521, 
3522, 3523, 3536, 3541, 3544, 3545, 3547, 3549, 3551, 3555, 3556, 3573, 
3575, 3578, 3581, 3582, 3587, 3589, 3594, 3595, 3599, 3604, 3609, 3611, 
3620, 3628, 3634, 3635, 3636, 3639, 3647, 3652, 3656, 3658, 3663, 3666, 
3668, 3680, 3685, 3686, 3687, 3689, 3691, 3694, 3695, 3696, 3706, 3709, 
3710, 3713, 3716, 3717, 3727, 3728, 3735, 3739, 3740, 3742, 3743, 3744, 
3751, 3764, 3772, 3775, 3776, 3778, 3782, 3787, 3789, 3792, 3794, 3798, 
3803, 3804, 3806, 3808, 3810, 3815, 3818, 3821, 3830, 3834, 3837, 3838, 
3842, 3847, 3849, 3850, 3863, 3867, 3868, 3883, 3884, 3889, 3890, 3899, 
3900, 3902, 3903, 3908, 3912, 3914, 3918, 3920, 3925, 3929, 3942, 3943, 
3944, 3947, 3949, 3956, 3957, 3960, 3962, 3963, 3978, 3980, 3981, 3986, 
3988, 3996, 4008, 4016, 4017, 4019, 4027, 4030, 4033, 4035, 4038, 4039, 
4044, 4049, 4057, 4072, 4073, 4075, 4077, 4082, 4086, 4092, 4106, 4111, 
4112, 4124, 4126, 4128, 4134, 4140, 4145, 4146, 4147, 4148, 4150, 4153, 
4157, 4162, 4176, 4177, 4179, 4188, 4189, 4190, 4193, 4195, 4203, 4206, 
4212, 4214, 4217, 4218, 4225, 4227, 4232, 4238, 4239, 4244, 4245, 4246, 
4247, 4257, 4262, 4271, 4273, 4279, 4283, 4288, 4291, 4293, 4300, 4308, 
4313, 4316, 4320, 4322, 4324, 4326, 4331, 4334, 4335, 4337, 4339, 4342, 
4346, 4349, 4350, 4351, 4352, 4354, 4357, 4361, 4376, 4377, 4386, 4388, 
4395, 4399, 4405, 4410, 4411, 4413, 4425, 4426, 4428, 4447, 4449, 4451, 
4453, 4459, 4471, 4475, 4480, 4488, 4493, 4496, 4501, 4504, 4506, 4512, 
4515, 4518, 4522, 4523, 4532, 4548, 4552, 4556, 4561, 4574, 4580, 4585, 
4587, 4595, 4598, 4602, 4605, 4610, 4614, 4619, 4621, 4626, 4628, 4632, 
4634, 4639, 4649, 4651, 4657, 4663, 4665, 4674, 4681, 4684, 4697, 4702, 
4703, 4705, 4710, 4711, 4714, 4725, 4728, 4736, 4738, 4740, 4751, 4754, 
4755, 4761, 4765, 4770, 4779, 4781, 4814, 4816, 4826, 4829, 4831, 4844, 
4850, 4873, 4877, 4880, 4885, 4890, 4892, 4896, 4899, 4900, 4902, 4903, 
4905, 4908, 4911, 4914, 4915, 4919, 4924, 4933, 4937, 4944, 4955, 4960, 
4967, 4969, 4974, 4992, 4994, 4997, 5000, 5002, 5003, 5007, 5019, 5021, 
5028, 5029, 5035, 5039, 5050, 5074, 5080, 5083, 5092, 5093, 5104, 5108, 
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5110, 5117, 5119, 5120, 5131, 5138, 5156, 5158, 5166, 5167, 5176, 5185, 
5192, 5194, 5200, 5203, 5219, 5225, 5226, 5233, 5236, 5239, 5251, 5259, 
5265, 5270, 5271, 5274, 5285, 5288, 5290, 5313, 5316, 5318, 5325, 5328, 
5336, 5337, 5346, 5364, 5389, 5405, 5408, 5415, 5432, 5439, 5440, 5448, 
5457, 5460, 5463, 5464, 5466, 5468, 5470, 5474, 5482, 5484, 5530, 5549, 
5567, 5589, 5591, 5598, 5612, 5627, 5632, 5664, 5676, 5680, 5718, 5719, 
5720, 5724, 5738, 5744, 5747, 5748, 5749, 5752, 5764, 5802, 5806, 5824, 
5841, 5859, 5880, 5885, 5892, 5896, 5899, 5927, 5928, 5947, 5956, 5957, 
5966, 5998, 6006, 6014, 6017, 6031, 6035, 6038, 6040, 6054, 6055, 6075, 
6081, 6084, 6089, 6092, 6113, 6145, 6166, 6169, 6180, 6192, 6195, 6221, 
6237, 6265, 6272, 6273, 6286, 6303, 6305, 6310, 6313, 6320, 6356, 6367, 
6456, 6469, 6510, 6515, 6520, 6521, 6524, 6541, 6561, 6609, 6615, 6642, 
6677, 6767, 6771, 6831, 6852, 6883, 6918, 6945, 6985, 7056, 7093, 7201, 
7379, 7387, 7461, 7560, 7588, 7685, 7764, 7765, 7979, 8308, 8460, 8899, 
9079 and 9923. But these results are not very useful since border cells 
are present. 


§ Result on 700 cells 


Box size: 10 

No compression 

Number of cells: 700 

Number of vertices: 4380 

Number of cells in frame: 341 

Time for counting stats: 2191.39 seconds 

Number of faces connected to the first vertice at infinity: 
Time for finding cell volumes: 3.1500 seconds 
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Ne Ny ; Nein tCPU, stat. tcpu,A| tcpu,v 
700 4,380 181 341 | 2191.39 sec. | 137.13 sec. | 3.15 sec. 
minny, max Ny, ee Try. ony, g;fv,| hyn. | nv. 
10 46 : 46.784 6.8399 24.169 | 23.208 
Nv, 5y(Mv.) M3 (nv,) M*(nv,) Kny,c 
24 5.4023 A 1.4240x10? | 6.5877x10° 3.0184 
MIN Ny cin MAaXNy cin i ee Ov cin Vogcin Vincin Ny cin 
10 fi 40.974 6.4011 25.464 24.668 
Ny cin Ci 5 i M3 (nv ,cin) M4 (no ,cin) Kn cin 
26 : f 1.0338x10? | 4.917310 2.9461 
MINNR, ,cin y sCin n Onwgcee One, cin Teg Nv icin TRARNs scin NR, cin 
8 26 : 10.244 3.2006 15.809 | 15.495 
PN, cin bu(MR, cin ) M? (nx, cin) M2 (nx, cin) MA (nx, cin) Knxy cin 
16 2.5550 10.214 12.922 | 3.073310? 2.9461 
minny.,¢ | Maxnx,,c Arx.,c On Re On, Rese Tig.Xerc Three | Nee 
8 27 A 12.439 3.5269 15.279 14.894 
15 2.7792 A 22.761 | 4.887110? 3.1673 
MINNY, cin | MAX NY, ,cin TR. scm Tr.) sein On(Re),cin Thg.Reycin | MARescin | MRe,cin 
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; 15.495 
M? (nx. cin i 
12.922 | 3.0733x10? 
TH(R,),c Fn(Np),c TNs. 
12.274 3.5034 14.857 
M3 (nx,,¢ M?(nyy,c 
21.05 4.6791 102 
MIN Nx, ,c MAX NX, cin TRS cin 
26 15.495 
ARs «ci Su (MR; cin 
2.5550 
min A, max A, Ah,c 
0.59209 1.0502 
A 
1.1935 x10 
min A,,, Ah, cin 
0.59209 : 2.0530 
Acin M*(Acin KA, 
1.3800 2.0572 2.9782 3.1031 1.0800 
min Afy,c max A fre Afre Apne Agne Ag, fr,c Ah, fr,c 
8.5015x10—8 | 1.1941x10—> | 4.8942x1077 |2.9253x1077 | 3.3091x10—> | 7.2835x1077 |1.5079x1077 








NX 5,c 


RR 5 ,cin 
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Afne 5u(Asre) M? (Afric) M3(Afrc) M* (Afr) KA gre 
1.7137x10-5 | 5.7986x10-7 |4.0635x10~-7 | 2.9478x10-7 |1.9816x10~7 | 2.9538x 10-7 
min Afrcin max A frcin Afnein TA yn a TA sr cin Ag, fr, cin Ah, fr, cin 
7.5700x107* |4.3579x 107° | 2.6048x 1073 | 1.3427x107% |5.1632x107% | 3.6182x10-* |2.6248x107% 
Afrycin Su (Afrein) | M?(Afrein) M?(Agrein) | M*(Asrcin) KA preci 
1.7644 107% |2.6301x10-% | 3.8077x10-° | 5.0023x10-% |3.9674x107° | 1.3808x 107 
min V, max V, Vv. oF, ov, Vqg,c Vhic 
0.88511 | 6.2052x10° | 2.7455x104 | 8.9303x10!° | 2.988410 22.548 6.3296 
Ve bu (Ve) M?(V,) M°(V.) M*(V,) Ky, 
7.6212 |_5.1126x10* | 8.917610" | 4.3202 x 10'” | 2.3781x10*4 |_ 2.990410? _ 
min V;,, oY. Ve, V9,Cin Vh, Cin 
0.88511 1.5264x104 | 1.2355x102 5.9195 4.4392 
Vein M3 (Vo.n) M*(Vo.n) KVen 
5.4355 3.9178x101° |_ 1.691510? 2 
min V¢y,¢ Vine Vg, fr,e Vh, frye 
4.6056 10-8 0.32288 | 1.4286x 10-3 1.5550x10-2 | 1.1733x10-° | 3.2935x 10-7 
Vere M? (Vor, M*(Vpr,c) EV prc 
3.9656x 10-7 | 2.6603x10~° |2.4145x107* | 6.0865x10-° |1.7433x10-> |_ 2.990410" ? 
min Verein max Vo rjcin Verein Frere Ve r,cin V9, fr, Cin Vh, fr, Cin 
1.2636x10~* 0.27064 | 2.9326x10—3 | 3.1108x10-* | 1.7638x10-7 | 8.4507x10—* | 6.3374x107* 


Vi rjcin 
7.7598x 10—4 





Su (Verein) 
3.9770x 10-3 


M?(Vir,cin) 
3.1017x 10-4 





M? (Vyr,cin ) 
6.6659x10—5 


M4 (Verein) 
1.6273 1075 





KV rein 


1.6915 x10 








Afrycin 


Ve 


Vorscin 
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§ AVS 


The first picture I created on AVS was a Trigonal Dipyramidal. The 
programme was the following .inp file. 


52000 

100 1.2910 

22 0 1.2910 

3 1 1.7321 1.2910 
4 10.5774 0 

5 1 0.5774 2.5820 
1itet i234 
21itet i235 
12 


The connection of modules was ReadUCD to ExternalEdges to 
UViewer3D. 


§ Problems related to ISD and NQS 


Some of the staffs at our Information Systems Department are un- 
professional in their jobs. Nicholas, for example, has once wiped out 
all about seven of my batch jobs on Network Queueing Systems. He 
emailed me on Monday 19" March that he had deleted one of my jobs 
because he thought it was doing nothing and consuming no CPU time. 
In his email he said it would not be able to reach its timeout limit in 
order to die out. 


I knew what he had written could be true and little minded the 
fact that he killed my task before bother telling me about it first. Two 
minutes after having read his email I found out at the ISD department 
that the number was not one but all, about seven in total, of my jobs 
had been deleted. Zaiem agreed at that point that it was not a right 
thing for a system administrator to do to delete users’ jobs without 
letting the owners know first. 


What happened in this case was this. Two weeks ago one of my 
NQS jobs correctly produced an input file of larger than 300 MB in size 
before it died due to lack of resource. I suddenly realised then that 
another one or two of my other jobs, if allowed to go on further, would 
produce output files of larger than 1 GB in size each. I tried to kill them 
using qdel. I found that I could not do so, so I asked Zaiem about it. 
For two weeks we tried everything in vain. We could not delete either 
one of them. Meanwhile, I did the only thing possible for me to do in 
order to prevent larger output files from them. Not being able to kill 
the job, I deleted their input files instead. 

It must have been for this reason that one of them went into a 
sleeping state as Nicholas said. But a sleeping programme could go on 


Vaen Sryayudhya, Editor January 2005 53 


Tyabandha Journal of Arts and Science Vol. 2, No. 1 


for a long time trying to read a nonexisting file without causing much 
trouble or even consuming much CPU time. 

When Nicholas came on the phone that morning, however, he com- 
pletely changed the story. Instead of saying that one of my jobs had 
been sleeping, it turned out to be that my jobs had messed up the whole 
NQS systems. He told me not to use NQS again after this. 

It is one thing to say that a job sleeps; it is another to say that it 
disrupts. In fact I think they are as opposite to one another as dull and 
vicious. 


I have just proved that two of my programmes which Nicholas 
killed actually completed very quickly when run interactively. I waited 
while the two ran in the foreground. One (f72.m) completed within 10 
minutes while the other (f71.m) took less than one hour to run. 


A programme submitted in January and killed in February due 
to its idleness had not gone to sleep as Nicholas suspected and Iain 
believed. The problem was with NOS. It has not been set up properly 
and it did not stop the process when after 24 hours. Instead, it let it 
run until well over 13 hours CPU time without doing anything about 
it. That process could, however, possibly have gradually slowed down 
for lack of resource, but not because diligence was missing. The people 
who looked after NQS obviously did not even realise this fact. Lack of 
time, as has been said, could not be a sound excuse since this was their 
major job. 

Some interesting facts I found out about NQS from this experience 
are the following. The command qdel can kill only jobs waiting in a 
queue. A running job needs to be killed from the platform it is running. 
NOS needs to be set up properly in order for it to kill a process which 
reaches its time limit. Matlab batch files might not expire when run 
on NQS. This queueing system is very convenient except that it needs 
absolute pathing in everything you do. Is there any way to configure 
NQS or write a shell script to help making it more user-friendly? 


§ More on number of vertices 3-d 

The following was done on a 3-d systems with a joggled-input 
option in qhull. The number of vertices, however, takes into account 
the boundary cells as well. The minimum n,,, was 6, maximum Ny, 
52, Tw,< 25.269, 2"¢-Moment 49.363, 3’¢-Moment 228.87, 4¢*-Moment 
8,704.9, on, 49.412, on. 7.0294, Ugn,,. 24.317, Mn,n,,. 23.365, Ty,¢ 24, 
6y(My,c) 5.5428, K(ny,c) 3.5725 . 

Following is a distribution graph. 
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Frequency versus number of vertices in a cell 
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Faces in different dimensions 


Considering only those cells bound withing the unit box, vertices 
and all, the following, namely Tables 6, 7 and 8, are the results from 
five simulations in two, three and four dimensions respectively, using 
the programme listed in Programmes section. 


n 


c 
tcpu(second) 


Random seed 
829247 134315 67453 432243 231215 


187 187 186 189 183 
170 167 168 163 170 
72 69 70 65 72 

169 165 166 159 170 


240 233 235 223 241 
20.36 20.39 20.09 20.6 19.76 


Table 6 Faces of Voronoi in two dimensions. N, = 100. 
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Distribution of number of vertices of nice cells 
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Distribution of number of vertices of nice cells 
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Distribution of number of vertices of nice cells 
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4 
Number 


(e) 


5 


6 7 


of vertices 


Random seed 


42398198 83250 


Ny 204 221 
Ny 148 147 

Ne 7 5 

ne 99 79 

pe 21.714 23.6 
(02). 4.5714 14.8 
m?(n?) 3.9184 11.84 
m?(n?) -4.6181 -16.128 
n2¢ 114 110 

fee 4.9211 4.7727 
(02) o¢ 1.2592 1.6268 
m*(n3, 1.2482 1.612 
m (nse 0.16453 0.99264 
ne 76 60 
(ue 5.0526 5.05 
(02)2! 1.1439 1.811 


m?((nb,)e) 1.1288 1.7808 
m3 ((n3,)c) 0.19004 0.28275 
nif 167 133 
tcpu(second) 24.28 25.67 
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34959 


224 
143 


743690 1321 


225 214 
155 146 

8 7 

109 97 

21 21.714 
34.286 21.905 
30 18.776 
171 69.831 
123 113 
4.8293 4.7788 
2.0772 1.656 


187 163 
34.23 26.46 
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Table 7 Various faces of Voronoi in three dimensions. N, = 50. 
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Distribution of number of vertices of bounded 2-d faces Distribution of number of vertices in 2-d faces of nice calls 























































































































Figure 36 Distribution of (ai) vc, (bi) vz¢, and (ci) v2! of the i** simulation 
on 3-d Voronoi. 


Random seed 
24098 802723 1453 732849 20480 


N, 1684 1719 1674 1586 1600 
Ny 938 921 860 889 904 
Ne 5 1 1 1 2 
ne 442 172 149 117 200 
pe 114 172 149 117 110 
(02). 1745 0 0 0 72 
2 (nv) 139.6 0 0 0 36 
3(n?) 547.2 0 0 0 0 
nsf 83 102 84 83 95 
je 8.8675 8.3137 8.5952 8.9398 8.7789 
(02) 3¢ 18.848 18.198 19.449 17.496 19.77 
m?(n¥,) 18.621 18.019 19.217 17.286 19.562 
m(n5,) 75.225 102.38 95.854 62.919 118.61 
nif 32 8 10 7 19 
( we 8.8675 9.25 10.8 11.143 9.8947 
(02) 18.848 13.643 21.511 14.476 28.655 
m?((n¥e)e) 18.621 11.938 19.36 12.408 27.147 
m “(in nee) 75.225 18.281 66.624 8.5364 195.1 
nf 11 0 5 3 7 
nif 2: 0 0 0 1 


tcpu(second) 315.07 358.94 324.56 281.85 283.82 


Table 8 Faces of Voronoi in four dimensions. N. = 100. 


The number of vertices in each 2-d face is a constant equals to 
three while that of a 1-d face is two. In the first run the number of 
vertices in each of the cells is 95, 108, 117, 120 and 130; in the last run, 
this number is 104 and 116. 
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Distribution of number of vertices of bounded 3-d faces 
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Distribution of number of vertices of bounded 3-d faces 
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Figure 37 Distribution of (ai) v3¢ and (bi) v2! of the i® simulation on 4-d 


Voronoi. 
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Table 9 contains the results obtained from a 4-d Voronoi network 
of 300 original nuclei. The numbers of vertices of the sixteen cells are 
97,99, 112, 141, 145, 160, 170, 171, 176, 176, 184, 186, 188, 192, 216 and 
235. 
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Figure 38 Distribution of (a) v3¢ and (bi) v2! of 4-d Voronoi, 300 nuclei. 


Table 10 is obtained from Voronoi in five dimensions. 
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Random seed 39378 Ny 16212 

Ny 6449 Nag 175 (o2)4# 352.2909 

Ne 1 be, 15.9200 m?((ne)e) 320.2645 

n® 864 (o2)a¢ 163.7752 m3 ((nip)e) 351.8362 

pe 864 2(n¥,) 162.8393 nif 1 

(02). 0 m3(ni,) 4.8764 x 10° (We 4 

m?(n®) 0 nit 11 a 1 

m3(n?) 0 (ue 29.0909 tcpu(second) 1.8908 x 104 


Table 10 Face statistics of 200 nuclei Voronoi in five dimensions. 
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Figure 39 Distribution of (a) vas and (bi) v* of 5-d Voronoi, 200 nuclei. 


Beam intersection study 


In this study of sectioning by a line the Voronoi in two dimensions, 
first generates on Matlab 500 points within a square box from —0.25 to 
1.25 in both axes. The beam is simply a straight line y = mz + c where 
m is the slope and ca constant. 
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Consecutive distances between consecutive intersections 
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Pencil y=3x-0.8 
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Figure 11 Intersection by a line. (a) is intersected by (b) y = 2a — 0.5; (c) 
is intersected by (d) y = —0.7z + 0.9; (e) and (f) are corresponding distances 
of respectively (b) and (d); (c) is intersected by (g) y = —5a + 4.9 and (h) 
y = 3a — 0.8. 


A natural basis for the normalisation is the expected distance. An- 
other possible basis is Wa = 0.064685. I call mean normalisation the 
normalisation using the first basis and homogeneity normalisation the sec- 
ond. Graphs of the closest pair distances look the same for both types 
of normalisation. 
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Simulation 1 2 3 4 
N. 500 1,000 
Ne 1,463 2,955 
Ne 239 463 
Ne 789 1,455 
Line equation y = 22-05 y 0.72+0.9  y 52+4.9 y=3x-0.8 
d 4.6867 x 10-2 3.163 x 10-7 3.2448 x 10-2 3.7843 x 107? 
a 3.9422 x 10-4 3.1174 x 1074 2.9005 x 10-4 4.9535 x 1074 
Normalisation 
by mean 120.4237 10.5582 10.5249 120.5881 
0.1713 0.3032 0.2656 0.3321 

—8.1332x 10-3 3.5783 x 10-2 —3.7762 x 10-3 2.2989 x 10-2 

by homogeneity 0.7245=£0.3070 0.68059=0.3799 0.6982=20.3665 0.81428=20.4789 
8.9936x 10-2 0.1404 0.1295 0.2202 


3.0936 x 10-3 1.128 x 10-2 —1.2853x 10-3 1.2412 x 10-2 
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A), where AB = B—A and CD = D—-C are vectors or directed lines and 

_ |CcD' AB' AB' AB' 
A, B,C, D are space vectors, r = CA! / Cp! CA! / Cp! 
Here AB = A+r(B-A), and CD = C+s(D-—C),0 <1r,s < 1 are directed 
lines. P exists ifO<r<land O<s<l. 


and s= 




















If the denominator is zero, then the two lines are parallel. 





B' 
CD' 





Also, if the nominator of r is zero, that is = 0, then both lines 





CD! 
CA’ 





are collinear. 

Consider the line section of Rayleigh distributed Voronoi where 
both the coordinates x and y are random numbers with Rayleigh distri- 
bution. The probability density function of the Rayleigh distribution is 

2 
y = f(a/b) = Bem. The mean of this distribution is b,/% and the vari- 
ance is 20°. With 6 = 1,2,...,1000 choose the random numbers from 
the Rayleigh distribution, then scale and use them as the coordinates. 

In Figure 41 the structure in (a) is intersected by (b) y = 3% — 0.8, 
(c) y=2, (d) y = 2z2, (e) y = 0.22, (f) y= —2 +1, (g) y= —@ 4+ 0.3. 
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Pencil y=3x-0.8, Rayleigh 
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x, Rayleigh 
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Pencil y=2x, Rayleigh 
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0.2x, Rayleigh 


Pencil y 














Vaen Sryayudhya, Editor 


January 2005 


82 


Vol. 2, No. 1 Tyabandha Journal of Arts and Science 


Pencil y=-x+1, Rayleigh 
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Pencil y=—-x+0.3, Rayleigh 
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Simulation 


Normalised 
by mean 


by homog. 


1 


1,000 

2,975 

988 

2,912 

y = 3a -0.8 
2.3621 x 1072 
8.1582 x 10-4 





121.2092 
1.4256 
5.5516 


0.7425. 0.8978 
0.7859 


2.2721 





yY=rr 
1.4253 x 107? 
6.2922 x 10-4 





121.7599 
3.0609 
28.351 


0.4480=L0.7885 
0.6144 


2.5493 





y = 22 
1.3853 x 1077 
2.9647 x 10-4 





121.2429 
1.5239 
7.6499 


0.4354=00.5412 
0.2890 
0.6316 





Table 12 Line intersection statistics of Rayleigh distributed Voronoi 


y = 0.22 
1.7434 x 107? 
3.0991 x 10-4 





121.0098 
1.0004 
1.6643 


0.548=00.5533 
0.3004 


0.2739 











5 

y 2+1 y= —-2£+0.3 

4.4401 x 10-2 1.027 x 10-7 

1.1584 x 10-3 3.5365 x 10-5 
120.7665 120.5791 





0.56204 
1.2114 


1.3956 1.0698 
1.0947 


3.293 








0.3245 
0.1755 

0.3228-£0.1869 
3.3813 x10—2 
5.9021 x10—3 
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The following codes generate and scale the points of Rayleigh dis- 
tribution. 


X=raylrnd([1:NumCell])’; 
Y=raylrnd([1:NumCel1])’; 
Max=0.8*max([X;Y]); 
X=X/Max; 

Y=Y/Max; 


The values of Cz and Dz will need to be adjusted manually from 
the pencil beam equation. It is the value of z at both points of inter- 
section between the beam and the [0,1] square box. A random Voronoi 
network can either be both homogeneous and isotropic or, nonhomo- 
geneous and nonisotropic depending on whether the probability distri- 
bution function is a constant. 


Percolation 


The 1970’s and the 1980’s saw a proliferation of variations on the 
theme of percolation. Every year there seemed to be a new percolation 
problem or two. For a Bethe lattice Chalupa et al (1979) reported a boot- 
strap percolation where those randomly occupied sites with less than m 
occupied neighbours are recursively emptied one by one until a stable 
configuration is reached. The problem they are interested in is that 
where the impurity concentration, dilution and crystal-field interaction 
compete in magnetic materials compete against the exchange interac- 
tion, resulting in the magnetic moments and consequently the magnetic 
order being destroyed. 

Wilkinson and Willemsen (1983) introduced invasion percolation. 
Working with Schlumberger they were interested in the real problem 
in the oil industry where water displaces oil via capillary action. Their 
approach was that of a constant flow rate, not one of a constant pres- 
sure as usual. Here water displaces oil in the smallest available pores. 
But when water completely surrounds any region of oil, no further ad- 
vance into that region will be possible, oil being incompressible. Such 
regions are called trappings and cause a problem generally known by 
the name of residual oil, an economic bane for oil industry. 

In the above example the hydrophobic versus hydrophilic prop- 
erty plays an important role in the replacement of oil in pores with 
water. And water is prevented from penetrating trapped oil regions by 
much the same principle as that which prevents the water in the con- 
tents of a sandwich from crossing the spread layer of butter to wet the 
hydrophilic bread. For this latter case the pores in question are those 
within the bread texture, and the soaking of the bread is best prevented 


86 January 2005 Vaen Sryayudhya, Editor 


Vol. 2, No. 1 Tyabandha Journal of Arts and Science 


provided that the trapped regions of butter or margarine percolates in 
two dimensions to form a single layer or film which entirely covers the 
sectioning surface of the bread. Moreover, there is a similarity also in 
the internal structure of both the rock and the bread. The density of 
pores in bread is not homogeneous as a result of the tension on the sur- 
face of the dough caused by internal air pressure originated from the 
yeasts inside, as well as because of the heat applied to it when inside 
the oven, which dehydrolises the surface and makes it dry and hard- 
ened. The same inhomogeneity can be found inside the rocks which 
form the oil reservoirs where the regions of oil are surrounded by rocks 
which are less porous and have less permeability. 


For Adler and Aharony (1988) a random walker, aka ant, treads on 
clusters. The ant enlarges a cluster by stepping on to an empty site next 
to it which meets certain conditions. They called this problem diffusion 
percolation. An example of a condition met is where the empty site has 
two or more occupied neighbours on a square lattice. 

Most percolation in studies happens by randomly toggling the 
phases of sites or bonds in regular lattices. Kerstein (1983) consid- 
ered randomly located spheres, take the complementary region of their 
union, and then perform percolation on the former. He showed that 
such percolation problem is equivalent to a percolation on a Voronoi 
lattice whose sites are the sphere centres. 

In the same way as an infinite loop in computer science means 
that one can not come to the termination of a programme going along 
the time dimension, an infinite cluster in percolation means that one 
can not come to the end of a cluster shifting along either one of the 
dimensions. The former case is along the one dimension of time and is 
only possible because of the time flows in one direction. Therefore half 
the line spanned is considered as being infinite. In one nondirectional 
dimension, except for the trivial case where one can consider the entire 
space as being one single cluster, no infinite cluster is possible. In 
simulation, when a cluster spans the whole of the space considered 
along any dimension we say that it is infinite since as one moves a 
cross section along that dimension, it always contains a section of the 
cluster. 

The bond percolation programme of two-dimensional Voronoi net- 
work by Tiyapan (1995, KNT4(iii), p. 78) takes O(n”) time provided we 
assume that the contribution made by the number of clusters together 
with the that by each cluster comparison amount to a linear time com- 
plexity term. In order to justify this assumption, consider first the 
number of clusters. The maximum number of clusters possible de- 
pends on the size of the system, in other words it must vary as n¥, 
where 0 < k&; < 1. This maximum value should be approximately 0.5 
because of symmetry between occupancy and void clusters. Next con- 
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sider the time involved in comparing two clusters. On Matlab this is a 
sparse vector comparison which is likely to involves some linked lists, 
and similarly should take time as a function of nk, 0 < ko < 1. Because 
on average the size of clusters is always small before Pc, kz will be 
less than 0.5. Therefore n*- nk < n and it is safe to assume O(n) time 
from both of them combined. q.e.d.. A C translation (Tiyapan 1995, ibid., 
p- 80) of this programme, though not as bad as it may seem because 
constant coefficients are small, gives O(n) time in comparison. 

In the field of geology Miller et al (2000) studies the analogy be- 
tween the dilatant slip in earthquakes and the hydrofraction occurred in 
melting and dehydration, the percolation of the latter in the permeabil- 
ity network internal to the fault being the cause of the former. Accord- 
ing to them, the existence of toggle switches in the permeability rules 
out the assumption that the permeability throughout the whole system 
is homogeneous. Instead, the system reorganises itself into sytems of 
different scales of interaction according to the degree and nature of its 
inhomogeneity. At the critical state the scale of interaction is equal to 
the scale of the model. 

The percolation probabilities given by Stauffer and Aharony (1994) 
are, the first number being for sites and the second for bonds, for the 
honey comb lattice 0.6962, 0.65271; square 0.592746, 0.50000; triangu- 
lar 0.50000, 0.34729; diamond 0.43, 0.388; simple cubic 0.3116, 0.2488; 
body-centred cubic 0.246, 0.1803; face-centred cubic 0.198, 0.119; 4-d 
hypercubic 0.197, 0.1601; 5-d hypercubic 0.141, 0.1182; 6-d hypercubic 
0.107, 0.0942; 7-d hypercubic 0.089, 0.0787. 


De Gennes et al (1959) investigated disordered binary solid solu- 
tion AB where active atoms A are randomly replace the nodes of the 
periodic matrix B. There exists a critical concentration A in B below 
which all clusters are finite, and above which both finite and non-finite, 
ie. infinite, clusters exist. Such solid solution in networks can repre- 
sent the spin waves in alloys with one ferromagnetic component or the 
impurity bands in semiconductors. They cited seminal work on perco- 
lation by Broadbent et al (1957), but no mention was made about the 
Ising model. 

In a way, percolation is similar to diffusion. In diffusion the parti- 
cles considered move about randomly, whereas in percolation they can 
only crop up randomly at predetermined locations on a network which 
is fixed. We could imagine, for instance, cars running along the roads 
within a traffic network as diffusing through them. Then the percola- 
tion could occur on a larger scale, that is the scale of a road. The cars 
move along, that means they diffuse; but the roads remain fixed, and so 
their phases could percolate. In other words, in diffusion the particles 
move while in percolation, whether there are moving particles or not, 
it is the phases that percolate. Since historically percolation began as 
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the study of diffusion of particles in a network of tubes in which the 
phases are naturally defined as the tubes being blocked or unblocked, 
these definitions have become most frequently used in other areas of 
application, for example in filtering membranes and traffic networks. 

But this is not necessarily the case. Instead of dealing with a fixed 
network, one may consider a model of percolation in a continuum, for 
example by randomly patching an area until all the patches connect 
with one another somehow and percolate. The patch could be of any 
shape, as well as polygonal and circular. We can consider the percola- 
tion in a certain area as having occurred when there appears a cluster of 
patches which traverses any two opposite sides. One application of this 
is in the study of occurrences of epidemics. Hoyle and Wickramasinghe 
(1979), having given a convincing argument in favour of viruses and 
various forms of diseases being carried to Earth from space by comets, 
talk about patchiness of pathogenic clouds. According to them, simul- 
taneous attacks across vast region rules out person-to-person theory. 
Moreover, influenza epidemics are generally characterised by sudden 
onsets and equally sudden ends. These epidemics and plagues may be 
thought of as the percolation of these patches in a sufficiently large, 
predefined area. The bacteria and viruses coming from space adding 
to gene give the possibility of jump patterns in evolution (cf Smith 
and Szathmary, 1995). The cells deliberately refuse to block viruses 
because they could prove to be useful in the long run, generally not 
by individuals but by the species. Historical examples are the disease 
described by Thucydides between 431 and 404 BC, five epidemics of 
‘English Sweats’ between 1485 and 1552, and more than ten influenza 
pandemic from 1700 to 1900. Example of diseases which are caused 
by bacteria and viruses are bubonic plague, chicken pox (varicella), 
cholera, common cold, Legionnaires’ disease, leprosy, measles, mumps, 
poliomyelitis, small pox, tuberculosis, and trachoma. Examples of ma- 
jor evolutionary transitions are those going from RNA to DNA, from 
prokaryotes to eukaryotes, and from asexual cloning to sexual propa- 
gation. Another example is the transition from primate to human both 
of whom differ from each other neurologically in the ability to use lan- 
guage and the power to conceptualise. One description of the Great 
Plague in London (Dickens, 1851). There was a rumour that a few peo- 
ple died in the winter of 1664. In May 1665 the disease burst out in St. 
Giles’s which raged through July and September in every part of Eng- 
land; approximately 10,000 people died in London alone. But then the 
equinox winds virtually blew the disease away, and the Plague quickly 
disappeared. The existence of interstellar organic matters is supported 
by strong evidences with more and more complex substances constantly 
found (cf Hoyle and Wickramasinghe, 1978). 

The growth mechanism of clusters in percolation is all so found 
in biology. Williams and Bjerknes (1972) simulate a tumour in the 
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basal layer of an epithelium. The basal cells become less sensitive to 
Charlone which controls cellular division, and thus they divide « times 
faster than the normal cells where & > 1 is the carcinogenic advantage. 
Abnormal cells interior to the basal layer divide, push, and then replace 
the neighbouring cells leaving the overall configuration unchanged ex- 
cept at the border where the abnormal cells exert a thrust of k — 1 on 
their normal neighbours, that is 4% = (« — 1)n where n and N are re- 
spectively the numbers of peripheral and total abnormal cells. They 
found that dimensionality of fractal is involved and the dimension 1.1, 
instead of 1, must be assigned to the periphery, which means that n 
is proportional to N°-°> not N°°. They found that abnormal cells push 
out faster in this order: the triangular, square, and hexagonal lattice. 
We may explain this, by looking at the coordination numbers of these 
three lattices, that the higher coordination number the lattice sites have, 
the slower the cluster expands. Added coordination means a higher 
degree-of-freedom the newly divided cells have to move about while 
still remaining local. 

There is no percolation in one dimension because in such case 
the percolating cluster must necessarily contain the whole space. But 
there are some applications where the critical blockage is important, 
for example the blockage of drainage grilles by pea shingle is one- 
dimensional. Here the critical amount of blockage depends on the crit- 
ical rate of flow which in turn depends on the amount of water and the 
rate of accumulation of water to be drained. The maintenance of gullies 
and grilles is done by cleaning, flushing and grit bucket emptying (cf 
Harrison and Trotman, 2002). 


Voronoi percolation, 2-d 


The percolation of Voronoi tessellation in two dimensions can be 
achieved by the following algorithm. 


generate random points; 
generate Voronoi tessellation and Delaunay triangulation; 
find vertices within the square box bounded by lower and upper 
bounds; 

find neighbours of all cells, bonds, vertices, and edges; 
for unit = cell, bond, vertice and edge do 

find number of units; 

permute list of units; 

clear cluster list 1, cluster list 2, set of resultant clusters; 

cluster percolated < false; 

for i = 1 to number of units do 

existing cluster joined = false; 
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for j =1 to number of clusters in cluster list 1 do 
if the j*" cluster contains the i‘ unit in permuted list do 


merge the i** unit into the j*" cluster ; 
existing cluster joined < true; 
end 
if existing cluster joined is true 
move the j*” cluster of cluster list 1 to cluster list 2 ; 
for k = 1 to number of clusters in cluster list 1 do 
if the k*” cluster in list 1 touches the cluster in list 2 do 
merge the former into the latter; 
else 
append the former to cluster list 2; 
end 
end 
test percolation of the cluster just updated; 
move cluster list 2 to cluster list 1; 
clear cluster list 2; 
break; 
end 
end 
if existing cluster joined is false 
create a new cluster of size one and append it to list 1; 
end 
append cluster list 1 to set of resultant clusters ; 
end 
end 


Algorithm 1 Network percolation in 2-d. 





Pc probabilities are found tobe f)1Pc = 0.511040.0856,  5"pp 
= 0.3095 + 0.0523, in Be = 0.7231 + 0.0616 and laa he = 0.6801 + 
0.0468. 

And coordination numbers are —.2436 7714.7. = 5.2320, — .2979 
eee, = 9 D022. 0259 Hien = 2.8617 and — .0382 hay ee = 3.8064 . 


8(2) 











Voronoi percolation, 3-d 


In three dimensions the procedure of finding p, is similar to that 
used in the 2-d case. Algorithm 2 is what I wrote and used in doing 
the simulation. I developed it to run on Matlab, and therefore do every 
thing in matrix form. As a consequence, there are various types of data 
all of which are matrices. These data structures can be summarised into 
the following. 
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A list is a one-dimensional matrix whose dimension is the number 
of its members. An index is a matrix whose members are either one 
or zero. These index matrices in two dimensions are like a map or 
an array. They are normally one or two dimensions, and map the 
relationship between the members of one set and another. A set is an 
m xn matrix every one of the members a;; of which is a matrix. A 
pair matrix is a square diagonal index matrix that maps among members 
within the same set. It contains the information about relationships 
like neighbourhood, group membership, connection, etc. A pair matrix 
or an index can also contain information other than the existence or 
membership flags, for example the edge length in the case of a vertices 
pair edge matrix. In a bond or an edge the mid point represents its 
position. 

The programme can be divided into three parts, creating and ar- 
ranging the Voronoi data, finding the neighbourhood matrix and the 
percolation simulation for p,. Putting the vertices of a face in order 
amounts to juggling from one end of each stick, ie. edge, to another. 
Edges are traced in one direction only until all the edges are success- 
fully linked head to tail. Two lists receive the result from the tracing, 
vertices which match go to one of them while those which do not is put 
in the other and recycled. 


Algorithm 2 Managing the Voronoi data in three dimensions. 


x + create random points; 
(va, Ca) + create Voronoi tessellation; 
t + create Delaunay tessellation; 
find v,, such that for all vg, 0 < vg < 1; 
find c, such that all vertices of cg are in vp; 
Ce Cp; 
for all do 
mp < find mid bond coordinates; 
b; + find bond length; 
endfor 
for all pairs of cells do 
fa < find shared vertices when either of the two cells is in cy; 
endfor 
for all vp, v < Ua} 
for all cn, f < fa; 
for all f do 
if it has three vertices then 
all the three possible pair combinations are neighbours; 
else 
(a,b,c) + find the face equation from three vertices; 
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6; < ka where k = 1/(a* +? +’) and i = 1, 2 and 3; 
max ¢ 0; 
for j = 1 to3 do 
if 9 < max then 
Im © Ji 
max + 4(jm)j 
endif 
endfor 
if jm =1 then 
POWIAHe 
elseif j,, = 2 then 
peu at 2; 
else 
pHUugry; 
endif 
d+ find delaunay triangulation from p and q; 
for all edges of all triangles of d do 
record the number of times they occur; 
endfor 
neighbours + vertices of edges which occur only once; 
for all f do 
order all their vertices; 
end 
endif 
endfor o 


The terms vertices and edges refer to the Voronoi tessellation only. 
The vertices and edges of the Delaunay tessellation are cells and bonds 
of the VT. The programme perco3d.m can find p, for all of these. For 
each one of them we must find the neighbourhood matrix as well as 
lists of the upper- and the lower boundaries. The order of blockages is 
completely decided in advance by finding a permuted list of numbered 
items. The programme finds the history of clusters for the whole range 
of p, i.e. from 0 to 1. 

To obtain the matrix of cell neighbours, first crop out between 5 
and 10 per cent the boundary in all directions. Then the neighbour 
matrix is found from the DT matrix. 

The bond neighbour matrix is simply the rearrangement of the 
cell neighbour matrix obtained. We have already found the vertices 
neighbour matrix earlier while searching for vertices shared between 
cells. This is rearranged for the use in the percolation programme, and 
then again for the edge neighbour matrix. 

Next is Algorithm 3 which carries out the percolation simulation. 
The cluster information is swapped alternately between three variables, 
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A, B and tmp, to facilitate the flow; A(i) is the i® cluster in A. All 
variables are assumed to be cleared at the start of the algorithm. At the 
end of each run the value of p, is computed, and the list p is reversed 
and reused for a second time if it was the first run. 


Algorithm 3 Voronoi percolation in three dimensions. 


for each permuted item p(i) in the list do 
joined + 0; 
for each cluster A(j) in A do 
if A(j) contains p(i) then 
A(j) — AG) N pli); 
B(1) + A(j); 
tmp < A; 
A+ imp — A(j); 
for all A(k) in A do 
if A(k) A B(1) then 
Bil) + A(k) N B(1); 
else 
B(+ +n) ¢ A(k); 
endif 
endfor 
percolated ?; 
A¢B; 
break; 
endif 
endfor 
if joined then 
A(+ +a) © pi); 


endif 
endfor a 
The resulted Pc probabilities from simulation are Bere: = 


0.2340+0.0448,  Y* .p, = 0.117840.0271, Y.p, = 0.2941+0.0831 
and Vp, = 0.4311 + 0.0324. 

And the coordination numbers are 3329 wwe = 11.1231 , 
5996 eas = 23.0548, .0245 Recs = 3.6926 and _.0432 Me a 
5.5002 . 














Percolation, 2-d Voronoi section 


The programme used does a 2-d section of the 3-d VT. Algorithm 4 
describes how it works. It assumes that the Voronoi tessellation already 
exists. Here d,, means a denominator while n, a numerator, V and C 
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means vertices and cells of the 3-d Voronoi tessellation while v and c 
those of the 2-d section. In particular, c € C 


Algorithm 4 Plane section of Voronoi in three dimensions 





(Ag, Ay, Az); — (@2 — %1, y2 — y1, 22 — 21)4 for all edges i; 
for all edges i do 
dm < ada + bAy + cAz; 
if d,, nonzero then 
Np € ax, + by + cz; 
t¢ —n-/dmi 
if 0<t<1 then 
(x,y, 2) + (a1 + Az, yi + Ay, 21 + Az)iz 
{us} + (a, y, 2); 
endif 
else 
atcat+eG 
dm < az + bAy + cAz; 
Np + ax, + by, + cz; 
t¢ —ny/dmi 
if 0<t<1 then 
(a, y, 2) & (41 + Aw, yr + Ay, 21 + Az); 
{us} & (x,y, 2); 
endif 
if n, =0 then 
{us} © (%1, 91; 21); 
{vs} < (2, yo; 22); 
endif 
endif 
endfor 
for all c, do 
find the Delaunay triangulation; 
count n, of all the triangles, add the numbers into a single list; 
endfor q 





The intersection of a plane ax + by + cz + d = 0 by the line which 
is defined by « = 2; + Act, y = y; + Ayt and z = z, + Azt, where Az, Ay 
and Az are respectively (x2 — 21), (yo — yi) and (z2 — 21), is determined 
by t = —(aa + by: + cz1 + d)/(aiAx + bAy + cAz). If the denominator is 
zero the line is parallel to the plane, and if the nominator is also zero 
contained therein. 

The results from the percolation simulation on the section V3 are 
VaQ3y = 05494401223, “#238 5, = 0.351540.0764,  Y823y, = 
0.7557 40.0757, xis" pe = 0.6210 + 0.0665. 
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The coordination numbers are 2212 pain ee = 4.6894 , .6021 
We ee 6.8005 0887 a ey 2 740by. ALB Aa, = 3,609! 3 

When sectioned by the plane (a,b, c,d) = (0.01, 0.5, 0.5, —0.5), our 
VT gives a picture as shown in Figure 42. Here the points shown are 
merely the average of the coordinates of all vertices in each cell. From 
ten simulations of these 118 cells, 300 bonds, 288 vertices and 405 edges 
having respectively the coordinate numbers of 5.0847, 9.7933, 2.8125 
and 3.7333, we obtain p, = 0.5220 + 0.0966, py, = 0.3107 + 0.0391, p, = 
0.7014 + 0.0573 and pe = 0.6259 + 0.0603. 

















When sectioned by the plane (a,b, c,d) = (0.5, —0.5, 0.5, 0.01), our 
VT gives a picture as shown in Figure 43. Here the points shown are 
merely the average of the coordinates of all vertices in each cell. From 
fourteen simulations of these 50 cells, 104 bonds, 140 vertices and 187 
edges having respectively the coordinate numbers of 4.1600, 8.2500, 
2.6714 and 3.5080, we obtain p, = 0.6914 + 0.1424, py = 0.4533 + 0.1108, 
Py = 0.7760 + 0.0821 and p, = 0.6929 + 0.0897.: 
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When sectioned by the plane (a,b, c,d) = (—0.7, —0.3, 1, 0.01), our 
VT gives a picture as shown in Figure 44. Here the points shown are 
merely the average of the coordinates of all vertices in each cell. From 
twenty simulations of these 121 cells, 305 bonds, 291 vertices and 411 
edges having respectively the coordinate numbers of 5.0413, 9.4426, 
2.8247 and 3.7518, we obtain p, = 0.5413 + 0.1107, py = 0.3584 + 0.0494, 
Py = 0.7077 + 0.0572 and pe = 0.6749 + 0.0470.: 




















For the purpose of these simulations, I have turned the code for 
finding percolation into a function. To my surprise and delight, this 
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same function works for both the 3-d VT and its 2-d section. In choos- 
ing a sectioning plane it is better if we choose the parameter d small, as 
the plane will then pass close to the origin. Also, choosing a+b+cx0 
seems to make a more wholesome section than otherwise. The codes for 
sectioning work well for oblique planes but do not like planes which 
are parallel to an axis. This shortcoming can be avoided if we make 
our plane only nearly parallel, when we want it to be parallel to an 
axis. Then to be able to view in a head on fashion such planes which 
have been plotted in three dimensions, we can look from the position 
(a,b,c), which is in effect the vector normal to the plane. The function 
mentioned above is listed as perc.m. 


Network percolation 


In the study of networks an important parameter is the coordina- 
tion number, which is the number of neighbours of an element which 
in the graph theory is usually the vertex. Each vertex or site of a graph 
is connected to each of its neighbouring vertices by a bond, so the co- 
ordination number of a graph is the number of bonds connected to a 
vertex. 


Clusters and their various characteristics play an important role in 
the study of percolation of networks. With a material science appli- 
cation in mind, Levy et al (1982) numerically represent the shape of a 
cluster by the shape parameter S, defined as S = b/N or more generally 
S = (1/2) 33_, i / 1 % where b or i is the number of bonds and N or 
v; the number of elements in a cluster having 8 or 7 bonds respectively. 


Programmes 


Network percolation, two dimensions 


VN=size(V,1); VCN=[]; Count=0; Xa=X; X=[]; 
Tmp=sparse(1,CNa) ; 
for i=1:CNa, 
Include=1; 
for j=1:VCNa(i,1), 
if (IXa(Ca{i}(1,j) ,1)==0) 


1% perco 

2 clear all; St=sum(100*clock); rand(’state’,St); CNa=200; 

3 Dim=2; X=rand(CNa,Dim); [Va,Ca]=voronoin(X); T=delaunayn(X) ; 
4 TN=size(T,1); VNa=size(Va,1); LB=0.05; UB=0.95; 

5 IXa=zeros(VNa,1); V=[]; Count=0; VCNa=(]; 

6 for i=1:CNa, 

7  VCNa=[VCNa;size(Ca{i},2)]; 

8 end 

9 for i=1:VNa, 

0 if ((Va(i,1)>LB & Va(i,1)<UB) & (Va(i,2)>LB & Va(i,2)<UB)) 
1 V=[V;Va(i,:)]; Count=Count+1; IXa(i,1)=Count; 

2 end 

3 end 

4 

5 

6 

7. 

8 

9 
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Include=0; 
end 
end 
if (Include==1) 
Count=Countt+1; C{Count,i}=[]; VCN=[VCN;VCNa(i,1)]; 
for j=1:VCNa(i,1), 
C{Count ,1} (1, j)=IXa(Ca{i}(1,j).1); 


end 
X=[X;Xa(i,:)]; Tmp(1,i)=Count; 
end 
end 
CN=size(C,1); T2=[]; T3=0); 
for i=1:TN, 
TmpA=(]; 
for j=1:3, 
if (Tmp(T(i,j))) 
TmpA=[TmpA, Tmp(T(i, j))1; 
end 
end 
TmpB=size (TmpA, 2) ; 
if (TmpB==2) 
T2=[T2;TmpA] ; 


elseif (TmpB==3) 
T3=[T3; TmpA] ; 
end 
end 
% for cells 
B=[]; BXX=sparse(CN,CN); NeCMat=sparse(CN,CN); Count=0; 
for i=1:size(T2,1), 
Count=Count+1; B=[B; [T2(i,1) ,T2(i,2)]]; 
BXX(T2(i,1) ,T2(i,2))=Count; 
BXX(T2(i,2),T2(i,1))=Count; NeCMat(T2(i,1) ,T2(i,2))=1; 
NeCMat (T2(i,2) ,T2(i,1))=1; 
end 
for i=1:size(T3,1), 
for j=1:Dim, 
for k=(j+1):(Dimt1), 
if (BXX(T3(i,j) ,T3(i,k) )==0) 
Count=Count+1; B=[B; [T3(i,j) ,T3(i,k)]]; 
BXX(T3(i,j),T3(i,k))=Count; BXX(T3(i,k) ,T3(i, j))=Count; 
NeCMat (T3(i,j) ,T3(i,k))=1; NeCMat(T3(i,k) ,T3(i,j))=1; 
end 
end 
end 
end 
BN=Count; A=X; N=size(A,1); LMat=sparse(1,N); 
UMat=sparse(1,N); LBc=0.2; UBc=1-LBc; 
for i=1:N, 
if (A(i,1)<=LBc) 
LMat (1,i)=1; 
elseif (A(i,1)>=UBc) 
UMat (1,i)=1; 
end 
end 
NeMat=NeCMat; Blocked=randperm(CN) ; 
% for bonds 
NeBMat=sparse (BN, BN) ; 
for i=1:CN, 
[a,b,c]=find(BXX(i,:)); nc=size(c,2); 
for j=1:(nc-1), 
for k=(j+1):nc, 
NeBMat (c(1,j) ,c(1,k))=1; NeBMat(c(1,k) ,c(1,j))=1; 
end 
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83 end 

84 end 

85 A=B; N=size(A,1); LMat=sparse(1,N); UMat=sparse(1,N) ; 
86 for i=1:N, 

87 if ((X(ACi,1),1)<=LBc) | (X(A(i,2) ,1)<=LBc)) 


88 LMat (1,i)=1; 

89 elseif ((X(A(i,1),1)>=UBc) | (X(A(i,2),1)>=UBc)) 
90 UMat (1,i)=1; 

91 end 

92 end 


93 NeMat=NeBMat; Blocked=randperm(BN) ; 
94 ¥for vertices 

95 Tmp=sparse(1,VN) ; 

96 for i=1:CN, 

97 for j=1:VCN(i,1), 

98 Tmp(1,C{i}(1,j))=1; 
99 end 

00 end 

01 Vv=O1; 

02 Count=0; 

03 for i=1:VN, 

04 if (Tmp(1,i)) 


05 Count=Count+1; Vv=([Vv;V(i,:)]; Tmp(1,i)=Count; 
06 end 
07 end 


08 VvN=size(Vv,1); 

09 for i=1:CN, 

10 for j=1:VCN(i,1), 

11 Cv{i}(1,j)=Tmp(1,C{it (1,9); 

12 end 

13 end 

14 E=[]; EVV=sparse(VN,VN); EVVv=sparse (VvN,VvN) ; 
15 NeVMat=sparse(VvN,VvN); Countv=0; Count=0; 
16 for i=1:CN, 

17. - Tmp=([Cv{i} (1,1: VCN(i,1)) ,Cv{i}(1,1)]; 

18 for j=1:VCN(i,1), 

19 Vi=Tmp(1,j); V2=Tmp(1, (j+1)); 


20 if (NeVMat (V1,V2) ==0) 

21 Countv=Countvti; NeVMat(V1,V2)=1; NeVMat(V2,V1)=1; 
22 end 

23 end 


24 «= Tmp=[C{i} (1,1: VCN(i,1)) ,c{i}(1,1)]; 
25 for j=1:VCN(i,1), 
26 Vi=Tmp(1,j); V2=Tmp(1, (j+1)); 


27 if (EVV (V1,V2)==0) 

28 Count=Count+1; E=(E;[V1,V2]]; 

29 EVV(V1,V2)=Count; EVV(V2,V1)=Count; 
30 end 

31 end 

32 end 


33 EN=Count; A=Vv; N=size(A,1); LMat=sparse(1,N); 
34 UMat=sparse(1,N); LBv=2*LB; UBv=(UB-LB) ; 

35 for i=1:N, 

36 «©6.if (ACi, 1) <=LBv) 


37 LMat (1,i)=1; 

38 = el seif (A(i,1)>=UBv) 
39 UMat (1,i)=1; 

40 end 

41 end 


42 NeMat=NeVMat; Blocked=randperm(VvN) ; 

43 {for edges 

44 NeEMat=sparse(EN,EN); MEV=sparse(EN, VN) ; 
45 [a,b,c]=find(EVV); nc=size(c,1); 
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46 for i=1:nc 


MEV(c(i),a(i))=1; MEV(c(i) ,b(i))=1; 


end 
for 


i=1:VN, 


a=find(MEV(:,1i)); 
if (~isempty (a) ) 


TmpN=size(a,1); Tmp=[a;a(1,1)]’; 
for j=1:TmpN, 
for k=(j+1) :TmpN, 
NeEMat (Tmp (1, j) , Tmp(1,k 
NeEMat (Tmp(1,k) ,Tmp(1,j 
end 
end 


))=1; 
))=1; 


end 


end 


A=E; N=size(A,1); LMat=sparse(1,N); UMat=sparse(1,N); 


for 


i=i:N, 


if ((V(A(i, 1) ,1)<=LBv) | (V(ACi,2) ,1)<=LBv)) 


elseif ((V(ACi,1),1)>=UBv) | (V(A(i,2) ,1)>=UBv)) 


LMat (1,i)=1; 
UMat (1,1)=1; 


end 


end 


NeMat=NeEMat; Blocked=randperm(EN) ; 
4% percol 
clear ClusA ClusB TSeries; NClusA=0; Perco=0; 


for 


i=1:N, 


Joined=0; 
for j=1:NClusA, 


if (ClusA{j,3}(1,Blocked(1,i)) ~=0) 


ClusA{j,1}=ClusA{j,1}+1; ClusA{j,2}(1,Blocked(1,i))=1; 
ClusA{j,3}=ClusA{j,3} | NeMat(Blocked(1,i),:); Joined=1; 


end 
if (Joined==1) 
for k=1:4, 
ClusB{1,k}=ClusA{j,k}; 
end 
NClusB=1; 
if (j==1) 
Tmp=ClusA; clear ClusA; 
for k=1:(NClusA-1), 
for 1=1:4, 
ClusA{k,1}=Tmp{(k+1) ,1}; 
end 
end 
elseif (j==NClusA) 
Tmp=ClusA; clear ClusA; 
for k=1:(NClusA-1), 
for 1=1:4, 
ClusA{k,1}=Tmp{k,1}; 
end 
end 
else 
Tmp=ClusA; clear ClusA; 
for k=1:(j-1), 
for 1=1:4, 
ClusA{k,1}=Tmp{k,1}; 
end 
end 
for k=j:(NClusA-1), 
for 1=1:4, 
ClusA{k,1}=Tmp{ (k+1) ,1}; 
end 
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209 end 

210 end 

211 for k=1:(NClusA-1), 

212 if (sum(ClusA{k,2} & ClusB{1,3}) ~= 0) 

213 ClusB{1,1}=ClusB{1,1}+ClusA{k, 1}; 

214 ClusB{1,2}=ClusB{1,2} | ClusA{k,2}; 

215 ClusB{1,3}=ClusB{1,3} | ClusA{k, 3}; 

216 ClusB{1,4}=ClusB{1,4} | ClusA{k,4}; 

217 else 

218 NClusB=NClusB+1; 

219 for 1=1:4, 

220 ClusB{NClusB,1}=ClusA{k,1}; 

221 end 

222 end 

223 end 

224 if ((sum(full(LMat & ClusB{1,2}))~=0) & ... 
225 (sum(full(UMat & ClusB{1,2}))~=0)) 

226 ClusB{1,4}=1; Perco=1; 

227 end 

228 NClusA=NClusB; ClusA=ClusB; clear ClusB; break; 
229 end 

230 end 

231 if (Joined==0) 

232 NClusA=NClusA+1; ClusA{NClusA,1}=1; 

233 ClusA{NClusA, 2}=sparse(1,Blocked(1,i),1,1,N); 
234 ClusA{NClusA, 3}=NeMat (Blocked(1,i),:); ClusA{NClusA,4}=0; 
235 end 

236 «© TSeries{i,i}=ClusA; TSeries{i,2}=Perco; 

237 end 


238 % Reverse 
239 Tmp=Blocked; Blocked=[] ; 
240 for i=1:N, 
1 Blocked=[Blocked, Tmp(1, (N-i+1))]; 
2 end 
3 % simulations 
4 Nc=0; TSnap=[] ; 
5 for i=1:N, 
246 if (TSeries{i,2}) 
? Nc=i; break; 
8 end 
9 end 
0 Pc=Nc/N; Cord=mean(sum(NeMat ,2)) ; 





Network percolation, three dimensions 


1% perco3d.m 

2 clear all; St=sum(100#clock); rand(’state’,St); 

3 CNa=300; Dim=3; X=rand(CNa,Dim); [Va,Ca]=voronoin(X) ; 
4 T=delaunayn(X); TN=size(T,1); VNa=size(Va,1); LB=0.05; 
5 UB=0.95; IXa=zeros(VNa,1); VCNa=[]; 

6 for i=1:CNa, 

7 VCNa=[VCNa;size(Caf{i},2)]; 

8 end 

9 MVCa=[]; 

10 for i=1:CNa, 

11 Tmp=ones(1,VCNa(i,1)); 

12 MVCa=[MVCa;sparse(Tmp,Ca{i},Tmp,1,VNa)]; 

13 end 

14 Vin=zeros(1,VNa); Count=0; 

15 for i=1:VNa, 

16 «if ((max(Va(i,:))<1) & (min(Va(i,:))>0)) 

17 Count=Countt+1; Vin(1,i)=1; IXa(i,1)=Count; 
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end 
end 
Tmp="Vin; Cin=ones(1,CNa) ; 
for i=1:CNa, 
if(sum(Tmp & MVCa(i,:))) 
Cin(1,i)=0; 
end 
end 
C=[]; count=0; VCN=[]; 
for i=1:CNa, 
if (Cin(i)) 
count=count+1i; TmpN=size(Ca{i},2); 
for j=1:TmpN, 
C{count ,1} (1,5) =IXa(Caf{i} (1, j)); 
end 
VCN (count ,1)=TmpN; 
end 
end 
CN=size(C,1); MidBCx=sparse(CNa,CNa); MidBCy=sparse(CNa,CNa) ; 
MidBCz=sparse(CNa,CNa); BLng=sparse(CNa,CNa) ; 
for i=1:TN, 
Tmp=(T(i,:),T(i,1)]; 
for j=1:(Dimt1), 
for k=(j+1):(Dimt1), 
if((Cin(1,Tmp(1,j)) | Cin(1,Tmp(1,k))) & ~BLng(j,k)) 
MidBCx (Tmp(1, 4) ,Tmp(1,4)) 2OC Ck, 1) #X(4,1))/25 
MidBCx (Tmp(1,k) , Tmp(1, j))=(X(k, 1)+X(j,1))/2; 
MidBCy (Tmp(1,j) , Tmp(1,k))=(X(k, 2)+X(j ,2))/2; 
MidBCy (Tmp(1,k) , Tmp(1, j) )=(X(k, 2)+X(j ,2))/2; 
MidBCz(Tmp(1, 3) ,Tmp(1,k))=(X(k,3)+X(j,3))/2; 
MidBCz (Tmp(1,k) , Tmp(1, j) )=(X(k, 3)+X(j,3))/2; 
dx=X(k,1)-X(j,1); dy=X(k,2)-X(j,2); dz=X(k,3)-X(j,3); 
TmpA=sqrt (dx*dx + dy*dy + dz*dz); 
BLng(Tmp (1,4) ,Tmp(1k))=TmpA; 
BLng(Tmp(1,k) , Tmp(1,j))=TmpA; 
end 
end 
end 
end 
Fa=[]; Count=0; FaC=[]; 
for i=1:CNa, 
if (Cin(1,i)) 
FaC{i,1}=0; FaC{i,2}=(]; 
end 
end 
for i=1:(CNa-1), 
for j=(it1):CNa, 
TmpA=0; TmpB=0; 
if (Cin(1,i)) 
TmpA=1 ; 
end 
if (Cin(1,j)) 
TmpB=1 ; 
end 
if(TmpA | TmpB) 
Tmp=MVCa(i,:) & MVCa(j,:); 
if (sum(Tmp) ) 
[a,b]=find(Tmp); Count=Count+1; 
Fa{Count ,1}=size(b,2); Fa{Count,2}=b; 
Fa{Count ,3}=[MidBCx (i,j) ,MidBCy (i, j) ,MidBCz(i,j)J]; 
if (TmpA) 
FaC{i,1i}=FaC{i,i} + 1; FaC{i,2}=[FaC{i,2},Count] ; 
FaC{i,3}{1,1}=i; FaC{i,3}{1,2}=j; 
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81 end 

82 if (TmpB) 

83 FaC{j,1}=FaC{j,1} + 1; FaC{j,2}=([FaC{j,2},Count]; 
84 FaC{j,3}{1,1}=i; FaC{j,3}{1,2}=j; 

85 end 

86 end 

87 end 

88 end 

89 end 


90 FaN=size(Fa,1); V=O; 
91 for i=1:VNa, 

92 if(Vin(1,i)) 

93 v=[V; [Va(i,:),i]]; 
94 end 

95 end 

96 VN=size(V,1); Tmp=sparse(VNa, 1); 
97 for i=1:VN, 

98 Tmp (V(i,4) ,1) =i; 

99 end 

00 F=Fa; 

01 for i=1:FaN, 


02. «for j=1: F{i, 1}, 

03 «FAG, 2} (1, j)=Tmp(Fa{i,2 (1,5) ,1); 
04 end 

05 end 

06 FN=FaN; FC=(]; 

07 count=0; 


08 for i=1:CNa, 
09 = if (Cin(i)) 


10 count=countt1i; FC{count,1}=FaC{i,1}; 

11 FC{count ,2}=FaC{i,2}; FC{count ,3}=FaC{i,3}; 
12 end 

13 end 


14 NghV=sparse(VN,VN); Tmp=F; TmpN=FN; 

15 for i=1:TmpN, 

16 TmpA=Tmp{i,2}; x=]; y=[]; z=[]; TmpB=Tmpfi,1}; 
17. «if (TmpB==3) 


18 for j=1:2, 

19 for k=(j+1):3, 

20 NghV(TmpA(1, j) ,TmpA(1,k))=1; 
21 NghV(TmpA(1,k) ,TmpA(1,j))=1; 
22 end 

23 end 

24 else 

25 for 1:TmpB, 

26 ile [V(TmpA (1,4), 1)]; 

27 ale V(TmpA(1,j), 2)]1; 

28 z=[z;V(TmpA(1,}) 31; 

29 end 


30 a=y (1) *(z(2)-z (3) ) +y (2) *(z(3) -2(1)) ty (3) (2 (1) -2(2)) ; 
31 b=z (1) * (x(2)-x (3) ) +z (2) * (x (3) -x (1) )+z(3) *(x(1)-x(2)); 
32 c=x (1) * Cy (2) -y (3) ) +x (2) * (y (3) -y (1) ) +x (3) # Cy (1) -y (2)) ; 





33 K=1/sqrt(a*a + b*b + ctc); 
34 Th{i}=K*a; Th{2}=K*b; Th{3}=K*c; Max=0; 
35 for j=1:3, 

36 if (Th{j}<Max) 

37 Max=Th{j}; jMax=j; 

38 end 

39 end 

40 if (jMax==1) 

41 Pry; 9=2; 

42 elseif (jMax==2) 

43 pHx; g=z; 
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t=delaunay (p,q) ; 
for j=1:size(t,1), 
for k=1:3, 
t(j,k)=TmpA(1,t(j,k)); 
end 
end 
Nt=size(t,1); TmpC=sparse(VN,VN) ; 
for j=1:Nt, 
TmpT=[t(j,:),t(j,1)]; 
for k=1:3, 
TmpD=sort ([TmpT(1,k) , TmpT(1, (k+1))]); 


TmpC (TmpD (1,1) ,TmpD (1,2) )=TmpC (TmpD (1,1) , TmpD (1,2) )+1; 


end 
end 
[k,1,m]=find(TmpC) ; 
for j=1:size(k,1), 
if (m(j)==1) 
NghV(k(j) ,1(j))=1; NghV(1(j) .k(j))=43 
end 
end 


end 


end 


Fed=[]; 


for 


i=1:FN, 


for j=1:2, 


Fed{i, j}=F{i, j}; 


end 


end 
for 


i=1:FN, 


Count=0; TmpN=Fed{i,1}; 
if (TmpN>3) 


Tmp=Fed{i,2}; TmpA=Tmp(1,1); Tmp=Tmp(1,2:TmpN) ; 
TmpN=TmpN-1; Count=Count+1; 
while (TmpN) 
a=TmpA(1,Count); TmpB=[]; Found=0; 
for j=1:TmpN, 
TmpC=Tmp (1,3); 
if(NghV(a,TmpC) & ~Found) 
TmpA=[TmpA,TmpC]; Found=1; 
else 
TmpB=([TmpB , TmpC] ; 
end 
end 
Tmp=TmpB; TmpN=TmpN-1; Count=Count+1; 
end 
Fed{i,2}=TmpA; 


end 


end 


LB2= 


2*LB; UB2=(UB-LB) ; 


% cells II 


Tmp= 


for 


ones(1,CNa) ; 
i=1:CNa, 


if ((max(X(i,:))>UB) | (min(X(i,:))<LB)) 


Tmp(1,i)=0; 


end 


end 


a=find(Tmp); TmpB=sparse(1,CNa); x=(; 


for 


i=1:size(a,2), 


TmpB(1,a(1,i))=i; x(i,:)=X(aCi),:); 


end 
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207 xn=size(x,1); TmpA=zeros(size(T)); TmpN=size(T,1); 
208 for i=1:TmpN, 
209 for j=1:4, 


210 «af ((Tmp(1,T(a, )))) 

211 TmpA (i, j)=TmpB(1,T(i,j)); 
212: end 

213 end 

214 end 


215 nghc=sparse(xn,xn) ; 

216 for i=1:TmpN, 

217. [a,b,c]=find(TmpA(i,:)); TmpB=size(c,2); 
218 if (TmpB>1) 


219 for j=1:(TmpB-1), 

220 for k=(j+1) :TmpB, 

221 nghe(c(1,j),c(1,k))=1; nghc(c(1,k) ,c(1,j))=1; 
222 end 

223 end 

224 end 

225 end 


226 A=x; N=size(A,1); LMat=sparse(1,N); UMat=sparse(1,N) ; 
227 for i=1:N, 
228 «=6if (A(i, 1) <=LB2) 


229 LMat (1,i)=1; 

230 elseif (A(i,1)>=UB2) 
231 UMat (1,1)=1; 

232 end 

233 end 


234 NeMat=nghc; Blocked=randperm(xn) ; 
235 % bonds II 
236 [a,b,c]=find(triu(nghc)) ; 
237 b=[a,b]; bn=size(b,1); Tmp=sparse(bn, xn) ; 
238 for i=1:bn, 
239 Tmp(i,b(i,1))=1; Tmp(i,b(i,2))=1; 
240 end 
1 nghb=sparse (bn, bn) ; 
2 for i=1:xn, 
3 a=find(Tmp(:,i)); 
4  if("isempty(a)) 
245 TmpN=size(a,1); 
6 for j=1:(TmpN-1), 
7 for k=(j+1) :TmpN, 
8 nghb(a(j,1),a(k,1))=1; nghb(a(k,1),a(j,1))=1; 
249 end 





250 end 
251 end 
252 end 


253 A=b; N=size(A,1); LMat=sparse(1,N); UMat=sparse(1,N); 
254 for i=1:N, 
255 if ((x(ACi,1),1)<=LB2) | (x(A(i,2) ,1)<=LB2)) 


256 LMat (1,i)=1; 

257. elseif ((x(A(i,1),1)>=UB2) | (x(A(i,2) ,1)>=UB2)) 
258 UMat (1,1)=1; 

259 end 

260 end 


261 NeMat=nghb; Blocked=randperm(N) ; 

262 4 vertices II 

263 [a,b,c]=find((NghV)); Tmp=sparse(1,VN); 
264 for i=1:size(a,i), 

265 Tmp(i,b(i,1))=1; 

266 end 

267 d=find(Tmp); TmpN=size(d,2); 

268 for i=1:TmpN, 

269 Tmp(1,d(1,i))=i; 
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270 end 

271 for i=1:size(a,1i), 

272, a(i,i)=Tmp(1,a€i,1)); b(i,1)=Tmp(1,b(i,1)); 
273 end 

274 nghv=sparse(a,b,c,TmpN,TmpN); TmpA=(]; 

275 for i=1:TmpN, 

276 TmpA (Tmp (1,d(i)),:)=V(i,1:3); 

277 end 

278 A=TmpA; N=size(A,1); LMat=sparse(1,N); UMat=sparse(1,N); 
279 for i=1:N, 

280 if (A(i,1)<=LB2) 


281 LMat (1,1)=1; 

282 elseif (A(i,1)>=UB2) 
283 UMat (1,i)=1; 

284 end 

285 end 


286 NeMat=nghv; Blocked=randperm(N) ; 
287 %4 edges II 

288 [a,b,c]=find(triu(NghV)); E=[a,b]; EN=size(E,1); 
289 Tmp=sparse (EN, VN) ; 

290 for i=1:EN, 

291 Tmp(i,E(i,1))=1; Tmp(i,E(i,2))=1; 
292 end 

293 NghE=sparse (EN ,EN) ; 

294 for i=1:VN, 

295 a=find(Tmp(:,i)); 

296 if (~isempty(a)) 


297 TmpN=size(a,1); 

298 for j=1:(TmpN-1), 

299 for k=(j+1) :TmpN, 

300 NghE(a(j,1),a(k,1))=1; NghE(a(k,1) ,a(j,1))=1; 
301 end 

302 end 

303 end 

304 end 


305 A=E; N=size(A,1); LMat=sparse(1,N); UMat=sparse(1,N); 
306 for i=1:N, 
307, if ((V(ACi, 1) ,1)<=LB2) | (V(ACi,2) ,1)<=LB2)) 


308 LMat (1,i)=1; 

309 elseif ((V(ACi,1),1)>=UB2) | (V(A(i,2),1)>=UB2)) 
310 UMat (1,i)=1; 

311 end 

312 end 


313 NeMat=NghE; Blocked=randperm(N) ; 

314 % percolation 

315 clear ClusA ClusB TSeries; NClusA=0; Perco=0; 
316 for i=1:N, 

317 Joined=0; 

318 for j=1:NClusA, 


319 if (ClusA{j ,3} (1,Blocked(1,i))~=0) 
320 ClusA{j,1}=ClusA{j,1}+1; 

321 ClusA{j,2}(1,Blocked(1,i))=1; 
322 ClusA{j,3}=ClusA{j,3} | NeMat(Blocked(1,i),:); 
323 Joined=1; 

324 end 

325 if (Joined==1) 

326 for k=1:4, 

327 ClusB{1,k}=ClusA{j,k}; 

328 end 

329 NClusB=1; 

330 if (j==1) 

331 Tmp=ClusA; clear ClusA; 

332 for k=1:(NClusA-1), 
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333 for 1=1:4, 
334 ClusA{k,1}=Tmp{(k+1) ,1}; 
335 end 
336 end 
337 elseif (j==NClusA) 
338 Tmp=ClusA; clear ClusA; 
339 for k=1:(NClusA-1), 
340 for 1=1:4, 
341 ClusA{k,1}=Tmp{k,1}; 
342 end 
343 end 
344 else 
345 Tmp=ClusA; clear ClusA; 
346 for k=1:(j-1), 
347 for 1=1:4, 
348 ClusA{k,1}=Tmp{k,1}; 
349 end 
350 end 
351 for k=j:(NClusA-1), 
352 for 1=1:4, 
353 ClusA{k,1}=Tmp{(k+1) ,1}; 
354 end 
355 end 
356 end 
357 for k=1:(NClusA-1), 
358 if (sum(ClusA{k,2} & ClusB{1,3}) ~= 0) 
359 ClusB{1,1}=ClusB{i,1}+ClusA{k, 1}; 
360 ClusB{1,2}=ClusB{1,2} | ClusA{k,2}; 
361 ClusB{1,3}=ClusB{1,3} | ClusA{k,3}; 
362 ClusB{1,4}=ClusB{1,4} | ClusA{k,4}; 
363 else 
364 NClusB=NClusB+1; 
365 for 1=1:4, 
366 ClusB{NClusB ,1}=ClusA{k,1}; 
367 end 
368 end 
369 end 
370 if ((sum(full(LMat & ClusB{1,2}))~=0) & ... 
371 (sum(full(UMat & ClusB{1,2}))~=0)) 
372 ClusB{1,4}=1; Perco=1; 
373 end 
374 NClusA=NClusB; ClusA=ClusB; clear ClusB; break; 
375 end 
376 end 
377. if (Joined==0) 
378 NClusA=NClusA+1; ClusA{NClusA,1}=1; 
379 ClusA{NClusA,2}=sparse(1,Blocked(1,i),1,1,N); 
380 ClusA{NClusA,3}=NeMat (Blocked(1,i),:); 
381 ClusA{NClusA,4}=0; 
382 end 
383 TSeries{i,i}=ClusA; TSeries{i,2}=Perco; 
384 end 
385 % Reverse 
386 Tmp=Blocked; Blocked=[] ; 
387 for i=1:N, 
388 Blocked=[Blocked,Tmp(1, (N-i+1))]; 
389 end 
390 % simulations 
391 Nc=0; 
392 for i=1:N, 
393 «if (TSeries{i,2}) 
394 Nc=i; break; 
395 end 
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396 end 
397 Pc=Nc/N; Cord=mean(sum(NeMat,2)); 


Network percolation, 2—d section 


%, section 
MVC=[]; 
for i=1:CN, 
Tmp=ones(1,VCN(i,1)); MVC=[MVC;sparse(Tmp,C{i},Tmp,1,VN)]; 
end 
CE=[(]; 
for i=1:EN, 
Tmp=MVC(: ,E(i,1)) & MVC(:,E(i,2)); 
if (sum(Tmp) ) 
TmpA=find(Tmp)’; TmpN=size(TmpA, 2) ; 
CE{i,1}=TmpN; CE{i,2}=TmpA; 
end 
end 
ie=[]; je=0; ke=(1; 
for i=1:EN, 
ie(i,1)=V(EGi, 2) ,1) -V(E(i,1) ,1); 
je(i,1)=V(E(i, 2) ,2)-V(E(i,1) ,2); 
ke(i,1)=V(E(i,2) ,3)-V(E(i,1) ,3); 
end 
20 a=1; b=.01; cc=0; d=-.5; v=[]; vC=[]; count=0; 
21 for i=1:EN, 
22. Tmp=(a*ie(i)+b*je(i)+cc*ke(i)); 


CANDUTUFPWNFOUWUANDUFLWNEH 





23 if (Tmp) 

24 vi=E(i,1); x1=V(v1,1); yi=V(v1,2); z1=V(v1,3); 
25 TmpA=(atx1tb*yl+cc*z1+d) ; 

26 v2=E(i,2); x2=V(v2,1); y2=V(v2,2); z2=V(v2,3); 
27 t=—-TmpA/Tmp; 

28 if ((t>=0) & (t<=1)) 

29 x=x1+(x2-x1)#t; y=yit+(y2-y1)*t; z=z1+(z2-z1)*t; 
30 count=count+1; 

31 v(count,:)=[x,y,z]; 

32 vC{count ,1}=CE{i,1}; vC{count,2}=CE{i,2}; 

33 end 

34 else 

35 if(“TmpA) % both nom and denom = 0 

36 count=countt1; 

37 v(count,:)=[x1,y1,z1]; vC{count,1}=CE{i,1}; 
38 vC{count ,2}=CE{i,2}; 

39 count=count+1; 

40 v(count,:)=[x2,y2,z2]; 

41 vC{count,1}=CE{i,1}; vC{count,2}=CE{i,2}; 

42 end 

43 end 

44 end 


45 vn=count; cC=sparse(CN,1); count=0; 
46 for i=1:vn, 

47 for j=i:vC{i,1}, 

48 if (~cC(vC{i,2} (j))) 


49 count=count+1; cC(vC{i,2}(j) ,1)=count; 
50 end 

51 end 

52 end 


53 cn=count; vc=vC; 

54 for i=1:vn, 

55 for j=i:vcfi,1}, 

56 ve{i,2}(1,j)=cC(vC{i, 2} (j)); 
57 end 
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for i=1:cn 


end 
for i=1i:vn, 
for j=l:vc{i,1}, 
c{vc{i, 2} (j) ,A}=c{ve{i,2}(j) ,1}+1; 
c{vc{i,2}(j) ,2}=[cf{vc{i,2}(j) ,2},i]; 
end 
end 


for i=1:cn, 


Tmp=[] ; 
for j=1:cfi,1}, 

tap: Erapa Catt 2 Cee tz 
en 
TmpA=min(Tmp,[],1); TmpB=max(Tmp, [],1); 
[TmpC , TmpD] =min (TmpB-TmpA) ; 
if (TmpD==1) 

TmpA=Tmp(:,2); TmpB=Tmp(: ,3); 


elseif (TmpD==2) 

TmpA=Tmp(:,1); TmpB=Tmp(: ,3); 
else 

TmpA=Tmp(:,1); TmpB=Tmp(: ,2) ; 
end 


TmpC=delaunay(TmpA,TmpB) ; TmpN=size(TmpC, 1) ; 
for j=1:TmpN, 
for k=1:3, 
TmpC (j,k) =Tmp(TmpC (j,k) ,4); 
end 
end 


TmpA=sparse(vn, vn) ; 
for j=1:TmpN, 
for k=1:2, 
for m=(kt1):3, 


Vol. 2, No. 


TmpA (TmpC (j,k) , TmpC (j ,m) )=TmpA (TmpC (j ,k) , TmpC (j ,m) ) +1; 
TmpA(TmpC (j ,m) , TmpC (j,k) )=TmpA(TmpC (j ,m) , TmpC (j,k) )+1; 


end 
end 
end 
[x,y,z]=find(TmpA); TmpB=[]; TmpC=[]; 
for j=1:size(x,1), 
if (2(j) ==1) 
TmpB=[TmpB;x(j) ,y(j)]; TmpC(y(j) ,1)=1; 
end 
end 
TmpA=(1; 
for j=1:size(TmpC,1), 
TmpA{j,i}=01; 
end 
for j=1:size(TmpB,1), 
TmpA{TmpB (j ,1) ,1}=[TmpA{TmpB (j,1),1},TmpB(j,2)]; 
end 
Tmp=Tmp(1,4); TmpB=Tmp; 
TmpC=sparse(Tmp,1,1,vn,1); count=c{i,1}-1; 
while (count>0) , 


if (~ (TmpC (TmpA{Tmp} (1) ,1))) 
Tmp=TmpA{Tmp}(1); TmpB=[TmpB,Tmp]; TmpC(Tmp,1)=1; 
else 
Tmp=TmpA{Tmp}(2); TmpB=[TmpB,Tmp]; TmpC (Tmp,1)=1; 
end 
count=count-1; 
end 
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40 








c{i,3}=TmpB; 
end 
for i=i:cn, 
Tmp=[0,0,0]; 
for j=i:cfi,1}, 
Tmp=Tmptv(c{i,2}(j),:); 
end 
Tmp=Tmp/c{i,1}; c{i,4}=Tmp; 
end 
Tmp=sqrt (a*atb*btcc*cc) ; 
de la/Tiip, b/ Tan, Gof taip) uzp=u; ux=[1,0,0]; Tmp=cross(u,ux) ; 
TmpA=sqrt (Tmp(1) *Tmp(1)+Tmp(2) #Tmp(2)+Tmp (3) *Tmp(3) ) ; 
uyp=Tmp/TmpA; uxp=cross(uyp,uzp) ; 
Relusp Ortyp Osae8, 050,050 11; 
vp=(R*[v’ ;ones(1,vn)])’; vp=vp(:,1:2); ad=min(vp,[],1); 
vp=vp- [ad(1) *ones(vn,1) ,ad(2) *ones(vn,1)]; caeL, 
for i=1i:cn, 
Tmp=R*[c{i,4}’;1]; Tmp=Tmp(1:2)’-ad; 
c{i,5}=Tmp; cs=[cs;Tmp] ; 
end 
LB=min(vp(:,1)); UB=max(vp(:,1)); Tmp=UB-LB; 
LBv=LB+0.1*Tmp; UBv=UB-LBv; Tmp=min(cs(:,1)); 
LB=min(cs(:,1)); UB=max(cs(:,1)); Tmp=UB-LB; 
LBc=LB+0.1*Tmp; UBc=UB-LBc; 
% cell 
cvm=sparse(cn,vn) ; 
for i=1:cn, 
for j=i1:cfi,1}, 
evm(i,c{i,2}(j))=1; 
end 
end 
nghc=sparse(cn,cn) ; 
for i=1:(cn-1), 
for j=(it1):cn, 
Tmp=find(cvm(i,:) & cvm(j,:)); 
if (“isempty (Tmp) ) 
TmpN=size(Tmp,2) ; 
if (TmpN>1) 
for k=1:TmpN, 
nghc(i,j)=1; nghc(j,i)=1; 
end 
end 
end 
end 
end 
N=cn; LMat=sparse(1,N); UMat=sparse(1,N) ; 
for i=i:cn, 
if (cs(i,1)<=LBc) 
LMat (1,i)=1; 


end 
if (cs (i, 1)>=UBc) 
UMat (1,i)=1; 
end 
end 
NeMat=nghc; Blocked=randperm(N) ; 
% bond 


[p.q,r]=find(triu(nghc)) ; 
b=[p,q]; bn=size(b,1); bcm=sparse(bn, bn) ; 
for i=1:bn, 
bem(i,b(i,1))=1; bem(i,b(i,2))=1; 
end 
nghb=sparse (bn, bn) ; 
for i=i:cn, 
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84 Tmp=find(bcem(: ,i)); 
85 if(~isempty(Tmp)) 


86 TmpN=size(Tmp, 1) ; 

87 for j=1:(TmpN-1), 

88 for k=(j+1):TmpN, 

89 nghb(Tmp(j) ,Tmp(k))=1; nghb(Tmp(k) , Tmp(j))=1; 
90 end 

91 end 

92 end 

93 end 


94 N=bn; LMat=sparse(1,N); UMat=sparse(1,N) ; 
95 for i=1:bn, 
96 if ((cs(b(i,1),1)<=LBc) | (cs(b(i,2) ,1)<=LBc)) 





97 LMat (1,i)=1; 

98 end 

99 if ((cs(b(i,1),1)>=UBc) | (cs(b(i,2) ,1)>=UBc)) 
200 UMat (1,1)=1; 

201 end 

202 end 


203 NeMat=nghb; Blocked=randperm(N) ; 

204 % vertice 

205 nghv=sparse(vn,vn) ; 

206 for i=i:cn, 

207. Tmp=[c{i,3},c{i,3}(1)1; 

208 for j=i:cfi,i}, 

209 nghv(Tmp(j) ,Tmp(j+1))=1; nghv(Tmp(j+1) ,Tmp(j))=1; 
210 end 

211 end 

212 LMat=sparse(1,vn); UMat=sparse(1,vn) ; 
213 for i=1:vn, 

214 «if (vp(i, 1) <=LBv) 


215 LMat (1,i)=1; 

216 end 

217. if (vp(i, 1) >=UBv) 

218 UMat (1,i)=1; 

219 end 

220 end 

221 N=vn; NeMat=nghv; Blocked=randperm(N) ; 
222 % edge 


223 [p,q,r]=find(triu(nghv)) ; 

224 e=[p,q]; en=size(e,1); evm=sparse(en,en) ; 
225 for i=1:en, 

226 evm(i,e(i,1))=1; evm(i,e(i,2))=1; 

227 end 

228 nghe=sparse(en,en) ; 

229 for i=i:vn, 

230 Tmp=find(evm(: ,i)); 

231 if (~isempty(Tmp)) 


232 TmpN=size(Tmp, 1) ; 

233 for j=1:(TmpN-1), 

234 for k=(j+1) :TmpN, 

235 nghe(Tmp(j) ,Tmp(k))=1; nghe(Tmp(k) ,Tmp(j))=1; 
236 end 

237 end 

238 end 

239 end 


240 N=en; LMat=sparse(1,N); UMat=sparse(1,N); 
241 for i=1:en, 
242 «if ((vp(e(i,1),1)<=LBc) | (vp(e(i,2) ,1)<=LBc)) 


243 LMat (1,i)=1; 

244 end 

245 if ((vp(e(i,1),1)>=UBc) | (vp(e(i,2) ,1)>=UBc)) 
246 UMat (1,1)=1; 
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247 end 
248 end 
249 NeMat=nghe; Blocked=randperm(N) ; 


Number of vertices 


1% numofvertices.m 

2 clear all; dimmin=2; dimmax=9; batches=5; 

3 dvn=[]; cpu=[]; nmax=1000; rand(’state’,sum(100*clock)) ; 
4 for i=dimmin:dimmax, 

5 for j=1:batches, 

6 n=round(nmax/i); x=rand(n,i); t=cputime; 

7 [v,c]=voronoin(x); cpu(i,j)=(cputime-t) /n; 

8 lend="floor(v); hend=~ (ceil (v)-ones(size(v))); 

9 lhend=lend & hend; in=min(lhend, [],2); 

0 dvn (i, j)=sum(in)/n; 

1 end 

2 end 

3 dvn=[(1:dimmax)’,dvn]; dvn=dvn(2:dimmax,:); 

4 figure(1); clf; 

5 for i=1:batches, 

6 semilogy(dvn(:,1),dvn(:,(it1)),’.’,’LineWidth’ ,2) ; 
7 hold on; 

8 end 

9 





edvn=[dvn(: ,1) ,sum(dvn(: ,2: (batchest+1)) ,2)/batches]; 

20 tmp=edvn(: ,2)./exp(edvn(:,1)); 

21 A=sum(tmp) /(dimmax-1); m=[dimmin,dimmax] ; 

22 semilogy(m, A*exp(m)) ; 

23 cpu=[(1:dimmax)’,cpu]; cpu=cpu(2:dimmax,:); figure(2); clf; 
24 for i=1:batches, 

25 semilogy(cpu(:,1),cpu(:,(i+1)),’.’,’LineWidth’,2); hold on; 
26 end 

27 ecpu=[cpu(: ,1) ,sum(cpu(: ,2: (batchest1)) ,2)/batches] ; 

28 tmp=ecpu(: ,2)./exp(ecpu(: ,1)); 

29 B=sum(tmp) /(dimmax-1); m=[dimmin,dimmax] ; 

30 semilogy(m, (B/35)*(exp(1)+2) .*m) ; 


Vertices per cell and cell ratio 


% numveachcell.m 
clear all; dimmin=2; dimmax=6; batches=5; nmax=3000; 
rand(’state’ ,sum(100*clock)) ; 
for i=dimmin:dimmax, 
for j=1:batches, 
n=round(nmax*2/i); x=rand(n,i); [v,c]=voronoin(x) ; 
fleet{i,j,i}=v; fleet{i,j,2}=c; fleet{i,j,3}=n; 
end 
end 
for i=dimmin:dimmax, 
for j=1:batches, 
v=fleet{i,j,1}; c=fleet{i,j,2}; n=fleet{i,j,3}; 
lend=~floor(v); hand (colt te) oncateice(sho): 
lhend=lend & hend; in=min(lhend, [] ,2); 
numve=[]; vcin=[]; 
for p=i:n, 
numvc=[numvc,size(c{p},2)]; flag=1; 
for q=1:numvc(p), 
if (~in(c{p}(q))) 
flag=0; break; 
end 
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22 end 

23 if (flag) 

24 vein=[vcin,numvc(p)] ; 

25 end 

26 end 

27 tmpn=size(vcin,2) ; 

28 rein(i,j)=tmpn/n; vec(i,j)=sum(vcin) /tmpn; 
29 end 

30 end 


31 dum=rcin; str={’c_{in} / c_fall}’}; 

32 dum=dum(2:dimmax,:); tmp=[]; 

33 for i=dimmin:dimmax, 

34 tmp=([tmp; i*ones (batches,1) ,dum((i-1),:)?]; 

35 end 

36 figure(1); 

37 semilogy(tmp(:,1),tmp(:,2),’.’,’LineWidth’ ,2); hold on; 

38 [p,s,mul=polyfit(tmp(:,1),tmp(:,2) ,4); 

39 x=(dimmin: .02: dimmax) ? ; y= polyval ( >x,[],mu); semilogy(x,y); 
40 dum=vec; dum=dum(2:dimmax,:); tmp= ah : 

41 for i=dimmin:dimmax, 

42 tmp= (tmp; i*ones (batches ,1) ,dum((i-1) ,dimmin: dimmax) ’] ; 

43 end 

44 figure (2) ; 

45 semilogy(tmp(:,1),tmp(:,2),’.’,’LineWidth’ ,2); hold on; 

46 edum=[dum(: ,1) ,sum(dum(: ,2: (batches+1)) ,2)/batches] ; 

47 tmp=edum(: ,2)./exp(edum(: ,1)); 

48 A=sum(tmp)/(dimmax-1); m=[dimmin,dimmax] ; 

49 semilogy(m, (A/70)*(exp(1)+4) .*m) ; 

50 elabeli Dimension”; Foutsiae? 143 ylabel(str,’FontSize’ ,14); 


Face statistics in n dimensions 


1% statsgenn.m by K N J Tiyapan, ist July, 2001 

2 echo off; clear all; format short g; more off; 

3 pti =fopen(’./v50.dat’,’r’); sci =fscanf(pti, ’%d’, 4); 
4 Dimension=sci(1,1); NumVA11 =sci(2,1); NumC =sci(3,1); 
5 sc2 =fscanf(pti, ’“%f’, (Dimension, NumVAl11]); 

6 VerticeAll =sc2’; CVMat =sparse(NumC, NumVA11) ; 

7 CFrame =ones(NumC, 1); VCFrame=zeros(NumVAl11,1); 

8 VFrame=ones (NumVA11,1); 

9 for i=1:Nunc, 

0 sci =fscanf(pti, *4d’, 1); 

1 for j=1:scl, 

2 sc2 =fscanf(pti, ’%d’, 1); 

3 Num =sc2+1; CVMat(i,Num) =1; 

4 if ( max(abs(VerticeAll(Num, :))) > 0.5 ) 

5 CFrame(i,1) =0; VFrame(Num,1)=0; 

6 end 

Z end 

8 end 





9 fclose(pt1) ; 

20 for i=1:NumC, 

21 VInC=find(CVMat(i,:)’); NumVInC=size(VInC,1); 
22 if (CFrame (i, 1)==1) 


23 for j=1:NumVInc, 

24 vCFrame (VInC(j,1),1)=1; 
25 end 

26 end 

27 end 


28 CVNiceCMat=[]; 
29 for i=1:NumC, 
30 if (CFrame (i, 1)==1) 
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CVNiceCMat=[CVNiceCMat ;CVMat(i,:)]; 
end 
end 
CNumVNiceCMat=sum(CVNiceCMat,2); NumV=sum(VCFrame) ; 
VVCFrameMat=zeros (NumV, 2) ; 
Vertice=zeros(NumV,Dimension); Count=0; 
for i=1:NumVAll, 
if (VCFrame (i, 1)==1) 
Count=Count+1; 
VVCFrameMat (Count ,1)=i; VVCFrameMat (Count ,2)=Count; 
Vertice (Count, :)=VerticeAll(i,:); 
end 
end 
pt3=fopen(’./n50.dat’,’w’); pt2=fopen(’./cb0.dat’,’r’); 
line=fget1 (pt2) ; 
scl =fscanf(pt2, ’4d’, 1); 
sc2=fscanf(pt2,’%f’, [Dimension,NumC]); Cell=sc2’; 
fclose(pt2); CNeighCCMat=sparse(NumC, NumC); t=cputime; 
FVA11Mat=(] ;sx FNumVAl1Mat-( . 
for i=1:(NumC-1), 
for j=(i+1) :NunC, 
VShared=and(CVMat(i,:), CVMat(j,:)); 
NumShared =sum(VShared, 2); 
NumFVA11Mat=size(FVA11Mat, 1); 
if (NumShared >= Dimension) 
CNeighCCMat (i,j) =1; CNeighCCMat(j,i) =1; Exist=0; 
for k=1:NumFVA11Mat, 


Mat chExistingFV=sum(and(VShared,FVA11Mat(k,:)),2); 


if (Mat chExistingFV>=Dimension) 
Exist=1; break; 
end 
end 
if (Exist==0) 
FVA11Mat=([FVA11Mat ; VShared] ; 
FNumVA11Mat=[FNumVA11Mat ;NumShared] ; 
end 


end 
FVMat=(]; FNumVMat=(]; FVCFMat=(]; FNumVCFMat=[] ; 
for i=1:NumFVA11Mat, 
VThisFace=find(FVA11Mat(i,:)’); 
NumVThisFace=size(VThisFace, 1); 
IncludeMe=1; IncludeMeToo=1; 
for j=1:NumVThisFace, 
if (VFrame (VThisFace(j,1) ,1)==0) 
IncludeMe=0; 
end 
if (VCFrame(VThisFace(j,1) ,1)==0) 
IncludeMeToo=0; 
end 
end 
if (IncludeMe==1) 
FVMat=[FVMat ;FVA11Mat(i,:)]; 
FNumVMat=[FNumVMat ; FNumVA11Mat(i,:)]; 
end 
if (IncludeMeToo==1) 
FVCFMat=([FVCFMat ; FVA11Mat(i,:)]; 
FNumVCFMat= [FNumVCFMat ; FNumVA11Mat (i,:)]; 
end 
end 
NumF VMat=size(FVMat,1); NumFVCFMat=size(FVCFMat ,1); 
FDim=Dimension-1; 
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fprintf(pt3,’Face dimension: %d\n’ ,FDim) ; 
fprintf(pt3,’Number of faces: %d\n’ ,NumFVMat) ; 
fprintf (pt3,’Number of vertices: \n [’); 
for i=1:NumFVMat, 
fprintf(pt3,’4d °’,FNumVMat(i,1)); 
if (mod (i, 10) ==0) 
fprintf (pt3,’...\n’); 
end 
end 
fprintf (pt3,’]\n’); 
fprintf(pt3,’Number of faces of nice cells: %d\n’ ,NumFVCFMat) ; 
fprintf(pt3,’Number of vertices: \n [’); 
for i=1:NumFVCFMat, 
fprintf(pt3,’%d ’,FNumVCFMat(i,1)); 
if (mod(i,10)==0) 
fprintf(pt3,’...\n’); 
end 
end 
fprintf (pt3,’]\n’); DVMat=FVCFMat; NumD=NumFVCFMat ; 
for d=3:Dimension, 
FaceCond=Dimension-d+2; DNeighDDMat=sparse(NumD, NumD) ; 
dVMat=[]; dNumVMat=[] ; 
for i=1:(NumD-1), 
for j=(i+1):NumD, 
VShared=and(DVMat(i,:), DVMat(j,:)); 
NumShared =sum(VShared, 2); NumdVMat=size(dVMat,1); 
if (NumShared >= FaceCond) 
DNeighDDMat (i,j) =1; DNeighDDMat(j,i) =1; Exist=0; 
for k=1:NumdVMat, 
MatchExistingdV=sum(and (VShared,dVMat (k,:)),2); 
if (Mat chExistingdV>=FaceCond) 
Exist=1; break; 
end 
end 
if (Exist==0) 
dVMat=[dVMat;VShared]; dNumVMat=(dNumVMat ; NumShared] ; 
end 
end 
end 
end 
FDim=Dimension-d+1; fprintf(pt3,’Face dimension: %d\n’,FDim) ; 
fprintf(pt3,’Number of faces: %d\n’ ,NumdVMat) ; 
if (FDim”=1) 
fprintf(pt3,’Number of vertices: \n [’); 
for i=1:NumdVMat, 
fprintf(pt3,’%d ’,dNumVMat(i,1)); 
if (mod (i, 10) ==0) 
fprintf(pt3,’...\n’); 


end 

end 

fprintf(pt3,’]\n’); 
end 
DVMat=dVMat; NumD=NumdVMat ; 
if (FDim==2) 

FVMat=DVMat ; 
end 

end 


Time=cputime-t; NumNiceC=sum(CFrame) ; NumVBound=sum(VFrame) ; 
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